% !Mode:: "TeX:UTF-8"
\documentclass[master,oneside]{seuthesis}
\usepackage{bm}
\usepackage{subfigure}
\usepackage{amsmath}
\usepackage{multirow}
\usepackage{slashbox}
\graphicspath{{D:/Resources/Papers/trunk/pic/}}
 % 这里是导言区

\begin{document}
\secretlevel{公开}    %学位论文密级分为"公开"、"内部"、"秘密"和"机密"四种
\studentid{091449}   %学号要完整，前面的零不能省略。
\title{迭代重建算法的研究与CUDA加速}{}{Study On Iterative Reconstruction and CUDA Acceleration}{}
\author{顾添锦}{Gu Tianjin}
\advisor{鲍旭东}{教授}{Bao Xudong}{Prof.}
                                % 可以不填
\degree{工学硕士} % 详细学位名称
\major[12em]{计算机科学与技术}
\defenddate{}
\authorizedate{}
\department{计算机科学与工程}{Department of Computer Science and Technology}
\duration{2009.09—2012.06}
\address{}
\maketitle

\begin{abstract}{锥束CT重建，迭代重建算法，SART，CUDA}

    锥束CT是当今国际CT研究领域最活跃的课题之一，与传统的二维扇束CT相比，它具有以下的几个优势：（1）具有更
    高的射线利用率，（2）扫描速度更快，（3）空间分辨率更高，这对工业无损检测以及医学影像诊断都具有重要的
    意义。

    锥束CT重建的算法主要分为解析重建法和迭代重建法，解析重建算法的经典近似算法是FDK算法，在大多数商用锥束
    CT中主要算法都是基于FDK重建算法。而迭代重建算法主要包括经典迭代重建与概率模型迭代重建，与解析重建相比
    迭代重建能够使用较少的投影数据重建较高质量的图像。但是迭代重建的速度较慢，需要叫多次迭代才能收敛。所
    以实际应用并不广泛。

    本文主要研究了三种经典迭代重建算法：代数重建算法（ART）、联合迭代重建法（SIRT）和联合代数重建法（SART
    ）。通过对三种算法的研究看到联合代数重建算法（SART）在重建质量和重建速度上都具有一定的优势，本文
    讨论了基于射线驱动和基于体素驱动的投影矩阵的计算方法以及各自对SART重建结果的影响，最后选择了基于射线
    驱动投影SART重建算法进行进一步的优化，文中把投影系统的旋转映射成重建区域数据的旋转，这样做避免每个角
    度投影矩阵的计算，进一步利用投影系统的对称性可以只需要计算1/4数量的射线所对应的投影矩阵，其余的3/4可
    以通过对称性计算所得，这样既减少了投影矩阵计算时间也减少了重建所需的内存，使得进行较高分辨率的重建成为
    可能。经过这些改进之后的SART算法在重建速度上有了一定的提高，但是算法上的改进之后的重建速度与解析重建法相比仍相去甚远。

    为了进一步提高SART的重建速度，本文中使用了NVIDIA公司推出的CUDA开发环境，有相对方便地将串行运行的程序该成在
    GPU上运行的并行程序，并且可以简单地通过硬件的性能的提升加快并行计算的速度，而无需为新的硬件修改代码。
    本文中基于CUDA加速的SART算法比CPU上运行的串行SART算法速度提高20至30倍，且图像质量与原图相比基本一致
    ，从而使得SART重建与传统的解析FDK算法重建有竞争的可能。
\end{abstract}
\begin{englishabstract}{Conebeam CT Reconstruction, Iterative Reconstruction,Simultaneous Algebraic
    Reconstruction Technology, CUDA}

    3D Cone-beam CT has been one of the most popular topics in the field of CT reconstruction research.It
    has several advantages compared with the traditional 2D fan-beam CT:(1)higher radiation
    effecient; (2) higher scan speed; (3)higher spacial resolution. These have important mean to
    industrial non-destructive testing and medical diagnosis.

    Algorithm of Cone-beam CT reconstruction mainly divide into two classes: iterative resontruction
    technology and analytical reconstruction technology. FDK algorithm is the most classic analytic
    reconstruction  technology, the algorithm of most commercial CT are base on FDK algorithm. Iterative
    reconstruction includes classic iterative reconstruction and iterative reconstruction based on
    probability model. Comparing with analytic reconstruction technology, iterative reconstruction can
    achieve the same reconstruction result with fewer projections. But iterativre reconstruction is slow,
    it needs iterate a few times to convergent, so it is not widely used.

    This article mainly researched on 3 classic iterative reconstruction technologies: Algebraic Reconstruction
    Technology(ART), Simultaneous Iterative Reconstruction Technology(SIRT) and Simultaneous Algebraic
    Reconstruction Technology. SART has advantage on reconstruction result and speed among the 3
    algorithms. This paper discussed two methods of computing the projection matrix: voxel-driven and
    ray-driven method. After compared these two methods, the ray-driven method was selected for further
    improvement. The rotation of the projection system can be mapped into the rotation of data of the
    reconstruction field, after this process there is no need to calculate the projection matrix at every
    angle. Further more in this stiuation it just needs calculate 1/4 of the ray in the projection matrix
    for the symmetry of the projection system. A lot of time and memory are saved by this operation, high
    resolution reconstruction becomes possible.

    For further improvement of the SART recontruction, this article presented an accelerated SART method
    by the enviroment of CUDA by NVIDIA Co,. With CUDA technology an sequential program can be easily
    change into an parallelized program run on the GPU, and the performance can be easily improved by
    changing better hardware without rewrite the code for new hardware. The CUDA accelerated SART can run
    20-30 times faster than the original method run on the CPU, and the reconstruction results are almost
    the same. This makes SART can  be comparable with the FDK series algorithm.
\end{englishabstract}
\begin{Main} % 开始正文
    \chapter{绪论}
    \section{课题背景及研究意义}
    计算机断层成像技术（Computed Tomography, CT）是通过对物体进行不同角度的射线投影而获得物体横截面信息的成
    像技术，1972年EMI公司的工程师G.N.  Hounsfied发布了世界上第一台商用CT扫描仪，使得在不进行介入手术的情
    况下可以对内脏进行检查\cite{mueller1998fast}成为可能。在这之后CT技术得到了长足的发展，直到今天他已经
    成为一种涉及发射物理学、数学、计算机、图形图像学和机械学等多个学科领域的高新技术,其具有诸多优点
    \cite{ymq2009}：(1) 非破坏性和非入侵性；(2) 非接触性；(3) 适用检测对象广泛，如生物体、各种金属和非金
    属材料；(4) 图像无影响重叠，密度分辨率和空间分辨率较高；(5) 可直接进行三维显示、图像处理、传输等数字
    处理。如今CT技术已广泛应用于诊断医学、无损检测和材料组织分析等工业领域\cite{ycsms}。

    CT扫描系统从问世到90年代初的20年间，随着计算机技术和性能的不断发展，CT系统不断朝着提高成像精度与质量
    的方向发展，但是主要集中于二维断层成像\cite{hj2006}。1998年多排螺旋CT系统多的出现被视为CT发展史上的一
    个里程碑。螺旋CT的射线源与多排探测器绕检测对象做螺旋扫描运动，多排探测器接受投影图像\cite{kalender20051}， 提高了扫描的效率。然而，目前常用的二维CT设备仍具有X射线利用率低，成像速度慢，
    难以对运动器官进行成像等缺点。所以锥束CT成为最近CT研究领域研究热点，锥束CT相较于二维扇束CT以及单排与
    多排螺旋具有以下一些优点\cite{Alex}：锥束CT具有更快的扫描速度和射线利用率，它在螺旋CT进行一次圆周扫描
    ，并有效减少X射线管的负载输出，降低扫描成本。除此之外，锥束CT可显著减少部分容积效应，能够获得比螺旋CT
    更高精度的空间分辨率，且其表现为X、Y、Z三个方向上各向同性。这对医学影像诊断和工业无损检车有着重要的意
    义。因此，锥束重建的研究成为目前CT重建的研究热点。

    随着这些CT技术的发展，研究人员提出了若干种重建技术，这些技术要被分成两大类\cite{mueller1998fast}：
    \begin{itemize}
        \item 解析法，基于傅里叶切片定理；
        \item 迭代法，把重建过程看作为解齐次线性方程组的问题。
    \end{itemize}
    相对于解析法，迭代重建算法具有很多优点，更适用于有限视角投影和投影图像噪声较为严重的情况，并且可以方
    便地根据具体成像条件引入与空间几何或测量值大小有关的约束条件，但是迭代重建需要较多的计算资源，所以在
    过去的商用CT领域一直未得到运用，解析重建算法是目前商业CT重建的主流算法\cite{wb2001cb}。但是随着计算机
    技术的不断发展特别的GPU通用计算技术的发展，经过并行化改进的迭代重建算法速度得到极大提高，使得迭代重建
    的临床应用成为可能。本文的主要工作就是在不影响迭代重建图像质量的情况下，利用CUDA对迭代重建算法进行GPU
    加速，提高迭代重建的速度，这对于迭代算法的工程应用以及CT重建具有较大影响。
    \section{国内外研究现状}
    \subsection{解析重建算法}
    解析重建算法的理论是依据1919年奥地利数学家Radon提出的Radon变换，根据Radon变换可以证明二维或三维物体可
    通过其无限投影重建起来。而对这个理论应用到诊断医学中的是Cormack，他在1963年的文章中不仅对此进行了可行
    性分析，并且还利用Fourier技术理论，给出了一种由无限多观测角度的投影值来求解目标切片的精确解析公式
    \cite{cormack1963}他也因此项研究与G.H. Hounsfield获得了1979年的诺贝尔生理学奖。这也说明了CT技术的出现
    对整个医学事业发展的影响。

    根据能否得到重建三维完备的Radon数据，锥束解析算法可以分为精确和近似两类算法。其中，根据重建物体区域不
    同，有可各自分为长物体重建和短物体重建。长物体重建算法是通过扫描物体的感兴趣区域（ROI）及附近的小区域
    ，仅仅利用纵向（Z轴）截断的锥束投影数据重建物体的ROI。相反，短物体算法就必须对整个物体进行扫描每次重
    建都必须对整个物体进行（图\ref{fig:1-1}）\cite{zk2003}
    \begin{figure}[htpb]
        \centering
        \subfigure[短物体问题]
        {
            \label{fig:subfig:1}
            \includegraphics[width=2in]{1-1-1.jpg} %文件名中不能出现括号
        }
        \hspace{0.5in}
        \subfigure[长物体问题]
        {
            \label{fig:subfig:2}
            \includegraphics[width=2in]{1-1-2.jpg} %文件名中不能出现括号
        }
        \caption{短物体和长物体重建示意图}
        \label{fig:1-1}
    \end{figure}
    在三维图像精确重建方面，1961年Kirillov就给出了锥束几何复值函数的逆变换公式，但复值锥束表示法不能在实
    际中应用。受Kirillov思想的启发，Tuy，Smith和Grangeat在1984年前后分辨提出了三种类型的锥束精确重建算法
    并给出了精确重建的条件——“每个与重建物体相交的平面，都和射线源扫描轨迹至少有一个交点”\cite{wb2001cb}。
    之后，研究人员在Tuy，Smith和Grangeat研究的基础上提出了更多高效的算法，其中最主要的是Grangeat类型的算
    法，此类算法易于在计算机上离散化，执行效率较高。1994年Defrise和Clack给出了一种采用位移可变的的滤波反
    投影精确重建算法\cite{defrise1994cone}。Kudo在1998年提出了真正通用的短物体精确重建算法，能利用截断的
    螺旋扫描投影数据进行精确的重建，成为了长物体精确重建算法的基础\cite{defrise1994cone}。2002年，
    Katsevich提出了螺旋轨道锥束精确重建\cite{katsevich2002analysis}。2004年，Zou和Pan通过变形Katsevich重
    建算法\cite{zou2004exact}，提出了反投影滤波（Backprojection Filtered，BPF）结构的重建算法。2005年，
    Pack等人通过引入Hilbert逆变换的计算算法，提出了基于R-line局部滤波的重建算法\cite{pack2005cone}。2007
    年，Katsevich对有正的曲率与挠率的光滑扫描轨迹和”圆+曲线“两类扫描轨迹情况下精确重建算法做出了重要贡献
    \cite{katsevich2007image}。国内上海交通大学庄天戈等提出了基于Radon变换的体积成像几何结构的完全性条件
    ，构造了基于中间函数的体积重建算法\cite{kanb2007con}。

    尽管精确重建算法能够重建出完美的结果，但是对采集环境、重建速度和数据的存储都需要较高的要求，所以至今
    尚无基于精确重建原理的商用CT出现。而这促使近似重建算法的发展与研究。近似锥束重建算法通常可以对不完整
    轨道投影进行重建、也可对探测角很少覆盖的对象进行重建、只需要一维滤波，正是这些优点，使得近似重建算法
    得到广泛的应用。

    在近似重建算法中，最著名的就是由Feldkamp，Davis和Kress在1984年提出的基于单圆轨道扫描轨迹的锥束重建算
    法，也就是FDK算法\cite{feldkamp1984practical}。之后，Grangeat、Defrise、Hu的研究成果表明FDK算法是精确
    重建算法在圆形轨迹上的特殊形式\cite{wy2008base}。但是FDK算法有其局限：只有当锥角较小时才能够去得较好
    的近似效果，并且只适用于圆形扫描轨迹，要求被重建的物体必须限制在球形视野区域内。1999年，Turbel提出
    P-FDK算法，将锥束数据重排成倾斜平行数据从而转化为二维重建问题\cite{turbell2001cone}，速度与FDK相当。
    2000年Grass等人基于P-FDK算法提出T-FDK算法，改善了FDK算法存在的低对比度衰减伪影。2001年，Grass又提出了
    HT-FDK算法，增加了可重建物体的体积。2004年，曾凯等人针对FDK算法中Radon平面数据缺失的问题，利用外推的
    思想提出了基于同心圆双轨道的CC-FDK算法\cite{zk2004base}，确保了在大锥角的情况下能够重建出较高质量的图
    像。Valton等人在2006年提出了适合于PET/CT成像系统的源点轨迹偏离中心的锥束重建算法\cite{valton2006fdk},
    该算法在偏离角小于17$^\circ$时比传统的FDK算法有更好的重建质量。
    \subsection{迭代重建算法}
    迭代重建算法又称直接重建算法，迭代算法中假设断层截面信息包含在一个离散的数据矩阵中，并且由若干的未知
    变量表示，然后由投影数据建立线性方程组，图像的重建过程在这里表现为方程组的求解过程。迭代重建算法主要
    分成两大类，一类是基于经典算法进行改进，另外一类是基于统计概率模型的迭代重建算法。

    经典的迭代算法包括ART（Algeriaic Reconstruction Technology）型以及SIRT（Simultaneous Iterative
    Reconstruction Technology）型。而ART又可以分为“加”型和“乘”等。Hounsfield的第一台CT实际上用的是ART算法
    \cite{ztg1992ct}。迭代重建的概念由Gordon、Bender、Herman等于于1970年首先引入图像重建领域。之后，
    Herman等人证明迭代重建算法是以“估值理论”作为其坚实的数学基础的\cite{herman1976iterative}。
    S.L.Wood(1987)又用随机滤波理论解释了迭代重建算法，并将迭代重建算法与变换法作了统一的处理
    \cite{slwood1978}。1984年Andersen和Kak对ART的体素更新方式进行了改进提出了SIRT，由逐个光线更新改为在所有光
    线都计算完成之后一次更新，从而有效地抑制条状伪影\cite{gilbert1972iterative}。但是SIRT算法的收敛速度比
    较慢，随后Andersen和Kak将ART算法与SIRT算法结合提出了Simultaneous ART算法（SART）
    \cite{andersen1984simultaneous}。

    第二类迭代重建算法从投影测量过程的随机性观点出发，把图像重建看成一个参数估计的问题，通过设计合理的目
    标函数而寻求是目标函数达到最优值的参数向量。八十年代初，Shepp\cite{shepp1982maximum}和
    Carson\cite{lange1984reconstruction}，在Poisson统计模型的基础上，分别提出了新的以最大似然（MLEM:
    Maximum Likelihood Expectation Maximized)为准则的重建算法，并由此开始了统计迭代重建的新领域，事实上，
    Shepp的重建算法是针对反射性断层（Emission Tomgraphy）投影数据而言的；而Carson则将其推广为穿透性断层（
    Transmission Tomography）投影数据的ML重建算法。受到最大似然的启发，后来有选取不同的先验概率模型利用最
    大后验概率（MAP:Maximum A Posteriori）重建图像\cite{nuyts1999simultaneous}。根据不同的概率模型，基于
    统计的迭代重建方法除了上面提到的两种还有有最小二乘法\cite{kaufman1993maximum}、最小范数法、Bayes估计
    、最小方差法和最大熵法等，采用优化准则随着优化理论的深入研究，已从最初单准则目标函数优化发展到多准则
    图像重建理论。
    \section{论文的主要工作和组织结构}
    迭代算法具有很多的优点，更适合于有限视角投影、噪声情况稍差等情况，本文讨论了迭代算法相对于解析算法的
    优势，并且研究了三种经典迭代算法各自的特点，最后选择SART算法作为本文的主要研究对象，并且利用Nvidia公
    司提供的CUDA开发环境对SART算法进行GPU加速算法的研究。其中主要的工作是对SART并行化算法的设计，经过改进
    后的SART算法课题利用有限的显存进行较高分辨率的重建。文章分为以下几个部分：

    第一章，绪论。对课题的研究现状的目的和意义作了综述，并且讨论了锥束重建目前的发展情况。

    第二章，CT重建的基础知识。主要介绍了CT成像的最基本原理，为后续算法的研究以及改进做理论准备。

    第三章，迭代重建算法的背景。从理论和试验的结果比较了三个经典迭代重建算法的各自特点，并且阐述为何选择
    SART算法作为本文的研究对象。

    第四章，SART算法的研究。主要讨论了SART重建算法的锥束实现过程。

    第五章，SART算法的CUDA加速。介绍了CUDA构架的硬件特点和程序执行模型，结合SART算法特点对SART算法进行改
    进实现并行加速。并给出结果分析。

    第六章，总结与展望。对本论文的研究内容进行总结，并对CT迭代重建算法的CUDA加速需要进一步解决的问题以及
    今后有待努力的方向进行展望。
    \chapter{CT成像的基本原理}
    CT图像是对设备采集的投影图像使用计算机按照一定的方案做数值计算而获得的。这种方案的目的是寻找投影数据
    与最终所需的重建物体密度之间的关系，并且根据这些关系得出重建区域的密度信息，即重建图像。而这里所说的
    重建方案就是各种CT重建方法。本章结合CT成像的物理和数学的基本原理，以及CT重建的重要理论基础——Fourier切
    片定理，二维CT的滤波反投影算法，以及迭代重建的思想来介绍CT重建算法的基本原理。
    \section{CT成像的物理基础}
    X射线是一种具有一定能量和穿透能力的高能电磁波，能够穿透一部分可见光不能穿过的物质（比如人体组织）。一
    般可见光的波长较长，光子在撞击到物体上的时候一部分发生反射，大部分被吸收；而X光的波长极短，光子所含能
    量较高，撞击物体之后，大部分能够穿越物体，其中部分光子因为光电效应、康普顿及瑞利散射使得X光强度发生衰
    减（图\ref{fig:2-1}），而衰减的规律服从Beer定理\cite{ztg1992ct}。X射线的穿透性与被照射物质等效原子序数
    和密度等信息相关，对于原子序数 较低组成的物质X射线的透射率较强，反之亦然。X射线的透射作用是CT成像的一
    个重要基础。
    \begin{figure}[htpb]
        \centering
        \includegraphics{2-1.jpg} %文件名中不能出现括号
        \caption{X射线透射物质示意图}
        \label{fig:2-1}
    \end{figure}

    考虑一束X射线射入均匀物质中，如图\ref{fig:2-2}所示。考虑Beer定理有：
    \begin{equation}
        I = I_0 \exp{(- \mu \Delta x) } \Rightarrow  \mu \Delta x = \ln{\left(\frac{I}{I_0}\right)}
        \label{eq:2-1}
    \end{equation}

    由式\ref{eq:2-1}可知，$\mu$值高的物体比$\mu$值低的物体使X光子衰减更多。例如，骨质的$\mu$ 比软组织的高，表
    明X光子穿透骨质比穿透软组织更加困难。另一方面，空气的$\mu$值几乎为0，说明在穿过空气的路径上，X射线输
    入和输出几乎没有改变。
    \begin{figure}[htpb]
        \centering
        \includegraphics[width=3in]{2-2.jpg} %文件名中不能出现括号
        \caption{X射线穿越均匀物质示意图}
        \label{fig:2-2}
    \end{figure}

    当X射线扫描的物质是不均匀的情况下，可将沿路径l分布的介质离散化成若干个连续的小块，当这些小块足够小时
    可以认为小块内部的介质是均匀的，具有相同的衰减系数，每个离散化小块的厚度假设为$\Delta x$，各个离散小
    块的衰减系数分别为$\mu_0$、$\mu_1$、$\mu_2$、$\mu_3$、\ldots、$\mu_n$，如图\ref{fig:2-3}。
    \begin{figure}[htpb]
        \centering
        \includegraphics[width=4in]{2-3.jpg} %文件名中不能出现括号
        \caption{X射线穿越不均匀物质示意图}
        \label{fig:2-3}
    \end{figure}

    X射线经过第一个小块之后的射线强度为$I_1$,经过第二块之后的射线强度为$I_2$ \ldots 最后出来的透射强度为$I_n$
    则
    \begin{align}
        I_1 = I_0 \exp{(- \mu_1 \Delta x) }  \label{eq:2-2} \\
        I_2 = I_1 \exp{(- \mu_2 \Delta x) } \label{eq:2-3}
    \end{align}

    将式\ref{eq:2-2}代入式\ref{eq:2-3}得：
    \begin{equation}
        I_2 = I_0 \exp{(- (\mu_1 + \mu_2)\Delta x) }
        \label{eq:2-4}
    \end{equation}

    沿着X射线传播方向继续累加衰减值，直到最后X射线离开被照射物体时的透射强度$I_n$：
    \begin{equation}
        I_n = I_0 \exp{(- (\mu_1 + \mu_2 +\ldots +\mu_n)\Delta x) }
        \label{eq:2-5}
    \end{equation}

    取式\ref{eq:2-5}的正值指数，并且表示为求和形式，得到
    \begin{equation}
        p = \sum\limits_{i=1}^{n}{\mu_i \Delta x }=\ln(\frac{I_0}{I_n})
        \label{eq:2-6}
    \end{equation}

    式\ref{eq:2-6}中的$p$即是投影。如果已知X射线的入射强度$I_0$和出射强度$I_n$，则根据式\ref{eq:2-6}得到
    一个以$\mu_i$（i=1,2,\ldots,n）为未知数的线性方程。当$\Delta x \rightarrow 0$时，式\ref{eq:2-6}
    可表示连续变化的求和，其积分形式为：
    \begin{equation}
        p =\ln(\frac{I_0}{I_n})= \int_{L}{\mu(l)}
        \label{eq:2-7}
    \end{equation}

    式\ref{eq:2-7}中，$\mu(l)$是衰减率关于路径$l$的连续函数。由投影求衰减系数函数$\mu(l)$的过程称为反投影。
    如果用二维密度函数$f(x,y)$来描述二维平面的衰减率，那么CT成像的问题就可以表述为：给定一个物体的被测量
    的线性积分，求各个点的衰减率，以产生二维密度数据。

    \section{解析重建算法}
    解析重建的理论基础是1917年由Radon提出的Radon变换。他从数学的角度证明可以用二维图像在各个角度的投影
    重建二维图像。本章主要介绍二维Radon变换和中心切片定理以简单阐述解析重建算法的理论依据。
    \subsection{二维Radon变换}
    二维Radon变换是一种直线积分的投影变换，它的定义形式有多种，本文以最常用的情况为例。假设在平面区域
    $R^2$中，任意一点$(x,y)$都可以用极坐标$(\rho,\theta)$表示，$\rho$表示直线到原点的距离，$\theta$表示
    直线$l$与$y$轴正向的夹角，函数$f(x,y)$ 为需重建图像。假设任意直线 $l\in R^2$，则函数$f(x,y)$的二维Radon变换定义为：
    \begin{align}
        &   p = \int_l f(x,y)dl=\int_l \hat{f}(r,\theta)dl  \notag   \\
        &  = \int_{-\infty}^{+\infty} \hat{f}(\sqrt{L^2+l^2} , \phi + tan^{-1}\frac{l}{L})dl
        \label{eq:2-8}
    \end{align}

    直线的方程在极坐标下可表示为$\rho = x \cos \beta +y \sin \beta$，则式\ref{eq:2-8}进一步可以表示为：
    \begin{equation}
        p =\int_{-\infty}^{+\infty}f(x,y)dl=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}f(x,y)
        \delta(x \cos \beta +y \sin \beta - \rho)dxdy
        \label{eq:2-9}
    \end{equation}
    其中，$(x,y)$表示重建像素点的在直角坐标系中的位置，$\delta(x)$表示采样函数。这个过程即是对图像在直线
    上的积分，根据之前X射线的物理原理，我们可以把X射线的衰减过程抽象成这个积分过程，X射线的投影值即是
    Radon变换在角度$\beta$的值。如果在多个角度进行投影，即可获得各个角度的Radon值，对这些获得的Radon值
    进行Radon逆变换，即可重建出原平面的二维图像$f(x,y)$。

    二维Radon逆变换的公式为：
    \begin{equation}
        f(x,y) = \hat{f}(r,\theta)=\frac{1}{2\pi^2}\int_{0}^{\pi}\int_{-\infty}^{+\infty}
        \frac{1}{r \cos{(\theta-\beta)}-\rho}\frac{\partial p}{\partial\rho} d\rho d\beta
        \label{eq:2-10}
    \end{equation}

    三维Radon变换是二维的推广，把二维中的线积分推广到面积分，每个面积分对应于Radon空间中的一个点。这个点
    就是该平面与该平面经过原点法线的交点。三维Radon变换空间就是由所有变换值组成。

    \subsection{中心切片定理}
    由于Radon逆变换公式在某些情形下是发散的，如在二维情形时，积分在$\rho=r\cos(\theta-\beta)$处发散，所
    以Radon逆变换不能直接用于反投影计算，需要用中心切片定理进行转换。中心切片定理是CT成像的重要理论基础。
    中心切片定理也被称为投影切片定理或者Fourier中心切片定理。它把投影的一维Fourier变换与图像的二维Fourier
    变换的一个切片联系起来。二维图像的中心切片定理证明：如果我们在不同的角度下获得足够多的投影数据，并
    对这些投影数据做Fourier变换，变换后的数据能够覆盖整个Fourier空间，一旦频率域函数$F(u,v)$取得所有值之后
    将其作一次Fourier反变换，就能得到重建图像$f(x,y)$。

    如图\ref{fig:2-5}所示即为Fourier重建方法，这也是常用的滤波反投影（FBP）重建算法的推导基础。
    \begin{figure}[htpb]
        \centering
        \includegraphics[width=5in]{2-5.jpg} %文件名中不能出现括号
        \caption{中心切片定理示意图}
        \label{fig:2-5}
    \end{figure}

    由图可知，对于每个投影数据先进行一维Fourier变换：
    \begin{equation}
        f(x,y) = \int_{-\infty}^{+\infty}P_{\beta}(t) \mathrm{e}^{-j2 \pi \omega t}dt
        \label{eq:2-11}
    \end{equation}

    式\ref{eq:2-11}中的$P_{\beta}(t)$是根据Radon变换计算$f(x,y)$在$\beta$方向上的投影：
    \begin{equation}
        P_{\beta}(t)=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}f(x,y)\delta(x \cos \beta +y \sin \beta - \rho)dxdy
        \label{eq:2-12}
    \end{equation}

    根据中心切片定理，可将式\ref{eq:2-11}的变换结果看成二维频率域中同样角度下过原点
    的直线上的值。在不同投影角度下所得的一维变换函数可在频域中构成完整的二维Fourier
    变换函数，将此二维变换函数做一次逆变换就得到所要求的空间域的密度函数$f(x,y)$，
    即重建图像。为了在二维逆变换中采用快速Fourier变换，通常在逆变换之前将极坐标形式
    的频域函数换成直角坐标形式的数据。

    于是我们就可以利用中心切片定理获得进行图像的重建的步骤：首先对各个角度上获得的投影数据进行Fourier
    变换，变换之后的数据将覆盖整个$(u,v)$平面。在函数$F(u,v)$获得全部值之后，对其进行一次二维Fourier
    反变换，这样就可以得到原始的密度函数$f(u,v)$，也就是重建图像。

    在Fourier重建的实际工程实现过程中，直角坐标系与极坐标之间插值转换时的精度问题是
    重建效果的关键。如果采用了不精确的插值方法，会导致伪影的存在从而限制了Fourier
    重建算法的应用场合。经过一系列的研究，科研人员尝试很多插值方法：最近邻点、双线性
    和截断sinc函数FIR插值器，并且证明了使用雅克比加权的二维周期sinc核函数的直接Fourier
    变换方法等价于FBP方法\cite{lh2008CT}。
    \section{迭代重建算法}
    滤波反投影（FBP）算法在实际运用中具有一定的局限性，比如它要求投影数据必须完全且均匀分布，滤波反投影
    的方程是连续形式的，实现时必须对图像进行离散化。在这种情况下，迭代算法是一种比较良好的选择。

    迭代重建算法的概念与解析重建算法最大区别在于前者把连续图像$f(r,\theta)$离散化。将整个图像区域划分
    为$J=n\times n$有限个像素，并以$\hat{f}(r,\theta)$表示。如图\ref{fig:2-6}表示射线经过离散化重建图像
    的过程：
    \begin{figure}[htpb]
        \centering
        \includegraphics[width=4in]{2-6.jpg} %文件名中不能出现括号
        \caption{离散化图像投影示意图}
        \label{fig:2-6}
    \end{figure}

    其中$x_1,x_2,\ldots,x_9$表示对应的像素值。由图可见各射线的和为：
    \begin{align}
        & p_1 = w_{11}x_1+w_{12}x_2+\dots+w_{19}x_9  \notag \\
        & p_2 = w_{21}x_1+w_{22}x_2+\dots+w_{29}x_9  \notag \\
        & p_3 = w_{31}x_1+w_{32}x_2+\dots+w_{39}x_9  \notag \\
        & \ldots                                     \notag \\
        & p_6 = w_{61}x_1+w_{62}x_2+\dots+w_{69}x_9
        \label{eq:2-13}
    \end{align}

    式\ref{eq:2-13}可以表达成更简洁的形式
    \begin{equation}
        p_i = \sum\limits_{j=1}^{9} w_{ij}x_{j} \quad i = 1,2,\dots,6
        \label{eq:2-14}
    \end{equation}

    或者矩阵表示形式
    \begin{equation}
        \bm{p} = \bm{Wx}
        \label{eq:2-15}
    \end{equation}

    式\ref{eq:2-15}中$\bm{p}=[p_1,p_2,\ldots,p_6]^T$,$\bm{x}=[x_1,x_2,\ldots,x_9]^T$，$W$为$6\times
    9$的矩阵。式\ref{eq:2-13}至\ref{eq:2-15}是根据9个像素，6条射线时候的特殊情况推得。我们再把他推广
    到一般情况，假设有$J$个像素
    $I$条射线的一般情况。此时：

    $\bm{x}$为$J$维矢量，称为图像矢量；

    $\bm{I}$为$I$维矢量，称为测量矢量；

    $\bm{W}$为$I \times J$维矩阵，称为投影矩阵。

    迭代重建的任务就是根据测量得到的$\bm{p}$和已知的投影矩阵$\bm{W}$（$\bm{W}$可以根据系统几何、焦斑形状、
    探测器响应以及CT系统其他物理参数确定）求：$\bm{x}$。\ref{eq:2-13}至\ref{eq:2-15}中的$w_{ij}$表示
    $i$号射线对$j$像素的加权因子，加权因子的计算有多种方法，最简单的情况根据射线经过是否经过像素点设
    加权因子为0或者1，如式\ref{eq:2-16}所示：
    \begin{equation}
        w_{ij} = \left\{
            \begin{aligned}
                & 1, \quad i\text{号射线经过}j\text{号像素；}  \\
                & 0. \quad \text{其他。}\\
            \end{aligned}
            \label{eq:2-16}
            \right.
        \end{equation}
        也可以把射线看成具有一定宽度，假设宽度为$\tau$（通常取像素宽度$\delta$）。这条粗线覆盖像素的一部分
        面积，覆盖面积所占像素的面积的比例就是该像素对所述射线投影的加权因子。例如$j$号像素的灰度值为
        $x_i$，$i$号射线与$j$号像素重叠区域面积为$\Delta s$，它与像素面积$\delta_2$之比为$w_{ij} = \Delta
        s / \delta_2$即为$i$号射线对$j$号像素的加权因子，具体投影矩阵的计算方法在之后的章节中有详细的阐述。

        有了以上基础之后，可以把CT图像的重建问题转化为线性方程组求解过程。最直观的办法是求矩阵$W$的逆矩阵
        $W^{-1}$，从而得到：
        \begin{equation}
            \bm{x} = \bm{W}^{-1}\bm{p}
            \label{eq:2-17}
        \end{equation}

        第二种解法是把所有经过$j$号像素的射线值累加，得到$j$号像素值：
        \begin{equation}
            x_i = \sum\limits_{i = 1}^{I}w_{ij}p_i, \quad j=1,2,\ldots,J.
            \label{eq:2-18}
        \end{equation}

        写成矩阵问题为：
        \begin{equation}
            \bm{x} = \bm{W}^T\bm{p}
            \label{eq:2-19}
        \end{equation}

        式\ref{eq:2-19}是反投影重建在离散像素情况下的形式，利用该式进行重建时伪影十分严重，但是式
        \ref{eq:2-19}有助于我们理解迭代重建算法。到目前为止，迭代重建遇到主要问题可能有：(1)一般像素数
        目及射线数目都极大直接求$\bm{W}^{-1}$难度较大，即使把$\bm{W}^{-1}$保存成稀疏矩阵，依然需要很大的
        计算量；(2)在一些情况中，投影个数远小于像素个数，线性方程组可能有无数组解；(3)在实际采集过程
        中，可能受物理偏差或者投影噪声等因素影响，无法得到重建结果。于是需要引入误差值，并估计一组解，使
        它在某一最优准则下达到最佳\cite{ztg1992ct}。于是可以对方程\ref{eq:2-15}做如下修正：
        \begin{equation}
            \bm{p} = \bm{Wx}+\bm{e}
            \label{eq:2-20}
        \end{equation}

        这里$\bm{e}$为误差矢量，可以为测量偏差和附加噪声，如探测器电子学噪声等。根据数值计算和最优化原理，
        估计的过程可由迭代实现，大体上包含如下步骤：（1）图像离散并且初始化；（2）选择迭代方式；（3）选取
        最优准则。迭代方式如上文所述，主要有经典的迭代（如ART、SIRT等等），基于统计的迭代（如EM、MAP等）。
        步骤三也有较多选择：最小二乘准则、最大均匀性及平滑准则、最大熵准则和贝叶斯准则等。
        \section{重建算法的评价标准}
        对于图像重建的算法主要包括两个部分，一是图像重建的质量，二是算法的可行性。而重建的质量主要包括空间
        分辨率，密度分辨率。而可行性主要指扫描和重建时间。
        \subsection{重建图像的评价标准}

        1. 主观评价

        图像质量的主观评价是由人的视觉感受对图像的质量进行评价。CT重建的目的就是显示被检测物体的断层图像以供人们进行
        观测，在实际医疗和工业应用中，投影重建图像的完全准确结果是未知的，没有参照对比主要以人的主观意识为评价
        标准，因此主观评价一直是评价CT重建图像质量的重要指标。但是主观评价自有其局限性，不同的检查个体对于同一
        图像的评价可能大相径庭，所以需要建立一套客观的评价标准以对重建图像的质量作客观真实的评价标准。

        2. 客观评价

        在采用模拟数据进行测试时，由于模拟数据的确切参数已知，所以可以将重建数据跟原始数据做精确的数值上的比较进而
        对重建质量作出客观评价。常用的模拟数据模型如Shepp-Logan标准头部模型，由许多大小密度不等的椭圆组合
        而成。常见的图像数值评价标准有：

        （1）图像相似性系数测量值$\varepsilon$
        \begin{equation}
            \varepsilon = \frac{\sum\limits_{i=1}^{N} (x_i-\bar{x}) (x_i^\ast -\bar{x^{\ast})}}
            {\left[ \sum\limits_{i=1}^{N}(x_i-\bar{x})^2 \sum\limits_{i=1}^{N}(x_i^\ast - \bar{x^\ast})^2 \right]^{1/2}}
            \label{eq:2-21}
        \end{equation}

        （2）归一化均方根距离测量值$d$
        \begin{equation}
            d=\left[\frac{\sum\limits_{i=1}^{N}(x_i-x_i^\ast)^2}{\sum\limits_{i=1}^{N}(x_i-\bar{x})^2}\right]^{1/2}
            \label{eq:2-22}
        \end{equation}

        （3）归一化平均绝对距离测量值$r$
        \begin{equation}
            r=\frac{\sum\limits_{i=1}^N |x_i-x_i^\ast|}{\sum\limits_{i=1}^{N}x_i}
            \label{eq:2-23}
        \end{equation}

        （4）图像信噪比$snr$
        \begin{equation}
            snr = 10 \times \lg\left(\frac{\sum\limits_{i=1}^{N}x_i^2}{\sum\limits_{i=1}^{N}(x_i-x_i')^2} \right )
            \label{eq:2-24}
        \end{equation}

        其中：

        $N$——重建图像中的像素个数；

        $x_i$——模型图像中第$x_i$号像素的灰度值；

        $x_i^*$——重建图像中第$x_i$号像素的灰度值；

        $\bar{x}$——模型图像中的平均灰度值；

        $\bar{x^*}$——重建图像中的平均灰度。

        以上的四种测量值强调了图像质量的不同方面。图像相似性系数$\varepsilon$反应的是重建图像与模拟
        图像之间的相似程度，$\varepsilon$越大说明两幅图像越相似，$\varepsilon$为1时表示两幅图像完全
        相同。归一化均方根距离测量值$d$较敏感地反应局部情况的误差，如果个别像素出现较大偏差会导致$d$
        较大。归一化平均绝对距离测量值$r$则较敏感地反应多数点小误差情况，与$d$恰好相反，它强调的是较多小误差
        的重要性，而不是少量大误差的重要性。信噪比$snr$衡量图像信号与噪声信号的比值，常常用分贝数表示。

        3. Profile曲线法

        Profile曲线法是一种更为直观的定量评价图像的方法，它直接比较基准图像和测试图像的对应像素的灰度值。
        具体做法是选择两幅图像中的同一行像素，在直角坐标系中按照图像的灰度变化，分别作出灰度变化曲线，
        比较两者之间的相似度，就可以对重建图像质量进行判断。
        \section{本章小结}
        本章主要介绍了CT重建的基本原理，从X射线的成像原理以及重建算法的理论基础。先介绍了X射线在射入均匀
        物质与非均匀物质时的数学模型，再通过这个数学模型讨论CT重建的方法。在解析重建算法中，介绍了解析重建
        算法的理论基础——Radon变换，以及Radon变换在实际运用中遇到的问题，从而引出中心切片定理以及它在实际
        重建中的运用——滤波反投影（FBP）算法。在迭代重建算法中，首先需要将图像离散化，根据投影环境的几何
        参数可得到迭代重建的投影矩阵，通过投影矩阵和投影图像可以建立一组线性方程组，迭代重建的过程就是
        对这个线性方程组求解的过程。最后本章给出了常见的图像重建评价方法，在之后的章节中将使用这些评价标准对各种重建算法
        的重建结果进行评价，进一步对各重建算法进行评价。
    \chapter{迭代重建算法的背景}
        在上一章中已经介绍了迭代重建算法的理论基础，本章分别介绍三种经典迭代重建算法，代数迭代重建算法（
        ART），联合迭代重建算法（SIRT）以及联合代数重建算法（SART）。本首先将介绍几种不同的投影矩阵计算方
        法，并且在介绍不同迭代重建算法时分别用基于不同模型的投影矩阵进行试验，比较各自投影矩阵对重建结果
        的影响。然后再分别介绍各个迭代重建算法的理论基础与具体实现方法，最后通过上章节中介绍的客观评价方
        法对各个重建算法的重建结果进行客观评价。
        \section{投影矩阵的计算}
        一个精确的投影矩阵对重建一副图像具有决定性的作用。在上节中，我们已经介绍过投影计算的基本情况
        ，本节以二维为例介绍投影矩阵具体计算方法。

        在所有计算投影矩阵方法中，最简单的模型是0,1模型，如图\ref{fig:3-1subfig:1}，当射线$p_j$经过像
        素$x_i$时设投影矩阵项$w_{ij}$为1，否则为0。图\ref{fig:3-1subfig:2}所示为基于长度加权的投影矩
        阵计算。投影矩阵项$w_{ij}$的值为像素$x_i$截取射线$p_j$长度的值。图\ref{fig:3-1subfig:3}表示基
        于面积加权的投影矩阵计算，这种模型中把射线当做具有一定宽度，投影矩阵项$w_{ij}$的值为像素$x_i$
        被射线$p_j$覆盖面积与像素面积的比值。
        \begin{figure}[htpb]
            \centering
            \subfigure[基于0,1的投影矩阵]
            {
                \label{fig:3-1subfig:1}
                \includegraphics[width=2in]{3-1-1.jpg} %文件名中不能出现括号
            }
            \subfigure[基于长度的投影矩阵]
            {
                \label{fig:3-1subfig:2}
                \includegraphics[width=2in]{3-1-2.jpg} %文件名中不能出现括号
            }
            \subfigure[基于面积的投影矩阵]
            {
                \label{fig:3-1subfig:3}
                \includegraphics[width=2in]{3-1-3.jpg} %文件名中不能出现括号
            }
            \caption{三种最基本的投影矩阵计算方法}
            \label{fig:3-1}
        \end{figure}
        在实际运用中，物体的内部衰减是连续的，在离散化的像素小格内部的衰减值实际上不是完全相等的。因
        此，重建出来的图像只是真实图像的离散化近似，重建出的图像可能会呈颗粒状。理论上我们可以通过使用插值
        方法重建出原始图像的连续描述。

        如图\ref{fig:3-2}所示，一条射线$p_i$射入重建区域中并且以等间距采样，插值核函数$h(u,v)$的中心
        放置在采样点上。再插值核函数范围中的所有重建像素都被累加，并且以插值核函数进行合适的加权。
        图\ref{fig:3-2}展示了在点$(X(s_{ik}),Y(s_{ik}))$处的采样值$s_{ik}$通过相邻的像素值计算。
        \begin{figure}[htpb]
            \centering
            \includegraphics[width=3in]{3-2.jpg} %文件名中不能出现括号
            \caption{在点$(X(s_{ik}),Y(s_{ik}))$处采样值$s_{ik}$插值计算示意图}
            \label{fig:3-2}
        \end{figure}
        $s_{ik}$的值通过式\ref{eq:3-1}计算：
        \begin{equation}
            s_{ik}=\sum\limits_{j}h(X(s_{ik})-X(x_j),Y(s_{ik})-Y(x_j))\times x_j
            \label{eq:3-1}
        \end{equation}

        射线$p_i$所对应的投影像素值$p_i$的值即为沿着射线所有的采样值$s_{ik}$累加：
        \begin{equation}
            p_i=\sum\limits_{k}\sum\limits_{j}h(X(s_{ik})-X(x_j),Y(s_{ik})-Y(x_j))\times x_j
            \label{eq:3-2}
        \end{equation}

        式\ref{eq:3-2}是式\ref{eq:3-3}的离散近似情况：
        \begin{equation}
            p_i=\int \left(\sum\limits_{j}h(X(s_{ik})-X(x_j),Y(s_{ik})-Y(x_j))\times x_j\right)ds_i
            \label{eq:3-3}
        \end{equation}

        将公式\ref{eq:3-3}进行重排得到：
        \begin{equation}
            p_i=\sum\limits_{j}x_j\int h(X(s_{ik})-X(x_j),Y(s_{ik})-Y(x_j))ds_i
            \label{eq:3-4}
        \end{equation}

        式\ref{eq:3-4}如图\ref{fig:3-3}所示，同式\ref{eq:2-14}类似，一个投影像素值$p_i$通过如下计算：
        \begin{equation}
            p_i=\sum\limits_{k}x_j\cdot w_{ij}
            \label{eq:3-5}
        \end{equation}

        因此权值通过沿着射线对插值核进行积分计算得到：
        \begin{equation}
            w_{ij}=\int h(X(s_{ik})-X(x_j),Y(s_{ik})-Y(x_j))ds_i
            \label{eq:3-6}
        \end{equation}
        \begin{figure}[htpb]
            \centering
            \includegraphics[width=3in]{3-3.jpg} %文件名中不能出现括号
            \caption{像素$x_{j}$与射线$p_{i}$的插值计算示意图}
            \label{fig:3-3}
        \end{figure}
        显然插值函数性能越好插值效果越好，加权值的计算也更加精确。图\ref{fig:3-4}展示了常用的插值核在空间
        域（图\ref{fig:3-4subfig:1}）和频率域（图\ref{fig:3-4subfig:2}）的波形。
        \begin{figure}[htpb]
            \centering
            \subfigure[各个插值核在空间域的波形]
            {
                \label{fig:3-4subfig:1}
                \includegraphics[width=2in]{3-4-1.jpg} %文件名中不能出现括号
            }
            \hspace{0.5in}
            \subfigure[各个插值核在频率域的波形]
            {
                \label{fig:3-4subfig:2}
                \includegraphics[width=2in]{3-4-2.jpg} %文件名中不能出现括号
            }
            \caption{各种插值核的波形}
            \label{fig:3-4}
        \end{figure}
        实验证明，线性插值核的效果已经非常接近于真实值，所以本章的实验中主要使用0,1值加权，长度加权已
        经线性插值核加权三种方法进行讨论。
        \section{代数重建法（ART)}
        代数重建算法（ART）是由Gordon，Bender和Herman为了解决三维物体重建的问题而提出的
        \cite{gordon1970algebraic}。这在当时滤波反投影算法应用还具有非常大局限性的情况下是非常巨大的突
        破。正如上文所说，Hounsfield发明的第一代CT扫描仪就是根据ART的原理制造的。尽管ART重建算法初始
        是为了重建三维物体。ART算法把重建的物体视为在平行投影结构下，一系列堆叠起来的二维切片的组合。
        在本章剩下来的讨论中，先把情况归约到二维的情况进行讨论，之后再扩展到三维的情况。

        ART可以写成线性几何问题：$\bm{Wx}=\bm{p}$，这里，$\bm{x}$是一个未知的（$N\times 1$维）列向量
        以保存在$n\times n \times n$大小的重建区域内全部$N=n^3$个体素。$\bm{p}$是一个$R$维列向量，
        $R$由每个投影的像素个数$R_m$乘以总共投影图像$P_\psi$集合中的投影个数$M$，即$R=R_m \times M$
        得到。$P_\psi$是指在一次扫描过程中所有投影图像的集合。$\bm{W}$是一个$R \times N$的投影矩阵，
        $\bm{W}$中的元素$w_{ij}$表示体素$v_{j}$对于射线$r_{i}$的影响。$\bm{Wx}=\bm{p}$可以写成于式\ref{eq:2-13}的线性
        方程组形式，正如上文所说直接求解这个线性方程组非常困难。所以我们这里介绍Kaczmarz求解这个线性
        方程组的方法。

        首先的问题是不能假设$R=N$，事实上在大多数应用中$R\ne N$的，有时候可以通过对投影图像使用插值算
        法强行使得$R=N$，但是随之带来的问题是投影矩阵$W$太过巨大无法保存。当$R>N$时最小二乘法可以用于
        解线性方程组，但是当$N$比较大时计算复杂性较大。在更一般的情况下，当$N>R$时方程$\bm{Wx}=\bm{p}$
        具有很多组解。所以ART重建算法的目标是通过投影图像获得重建区域最接近的估计值。

        Kaczmarz早在1937年就提出了解决这个问题的迭代方法\cite{kaczmarz1937}。在这个方法中，首先选择一
        个初始的体素向量，$\bm{x}=\bm{x^{(0)}}$，然后根据$\bm{x^{(0)}}$求一次近似图像$\bm{x^{(1)}}$，
        再由$\bm{x^{(2)}}$求二次近似图像，如此下去。ART校正方法被称为逐线校正，在根据$\bm{x^{(k)}}$求
        $\bm{x^{(k+1)}}$时需要求一个校正值$\bm{\Delta x^{k}}$，而$\bm{\Delta x^{k}}$只考虑一条射线投
        影的影响，只修正被射线$i$经过的那些体素。再修正完$i$号射线经过的体素之后，再处理$i+1$号射线，
        如此直到达到优化停止准则为止。因此，我们可以把每个体素的更新分解为三个阶段：(1)投影阶段，(2)
        矫正值计算阶段，(3)反投影阶段。

        重建区域$x$的更新过程的表达式如\ref{eq:3-7}所示：
        \begin{equation}
            x_j^{k+1}=x_j^k + \lambda \frac{p_i - \sum\limits_{n=1}^{N}w_{in}x_n^k}{\sum\limits_{n=1}^N w_{in}^2}w_{ij}
            \label{eq:3-7}
        \end{equation}

        式\ref{eq:3-7}中$\lambda$是一个松弛因子，他的取值范围在(0,1]区间内。但是在通常情况下，如果
        $\lambda$值过于接近于0的话会导致过度优化。将算法对\ref{eq:2-18}中的式子依次进行计算，完成之后
        可能某些重建像素不一定满足收敛条件，可以同样的方式进行下一次迭代。图\ref{fig:3-5}展示的是Kacmarz
        方法的几何过程，图中的两条直线可以表示成两个线性方程，图中所示个过程即为求解过程：先假设初始
        值$X^{(0)}$，将其投影到直线$L_1$上，接着将结果$X^{(1)}$投影到直线$L_2$上，之后再投影到
        $L_1$上，如此计算下去知道满足收敛条件，即可求得近似解。
        \begin{figure}[htpb]
            \centering
            \includegraphics[width=3in]{3-5.jpg} %文件名中不能出现括号
            \caption{Karmarz发求解过程示意图}
            \label{fig:3-5}
        \end{figure}
        图\ref{fig:3-6}是ART使用不同的加权模型在$z=100$处的重建结果，重建的参数为：$\lambda=0.25$，投影角度个数
        为90，迭代次数5次。
        \begin{figure}[htpb]
            \centering
            \subfigure[使用0,1加权]
            {
                \label{fig:3-6subfig:1}
                \includegraphics[width=2in]{3-6-1.jpg} %文件名中不能出现括号
            }
            \subfigure[使用长度加权]
            {
                \label{fig:3-6subfig:2}
                \includegraphics[width=2in]{3-6-2.jpg} %文件名中不能出现括号
            }
            \subfigure[线性插值核加权]
            {
                \label{fig:3-6subfig:3}
                \includegraphics[width=2in]{3-6-3.jpg} %文件名中不能出现括号
            }
            \caption{ART重建结果}
            \label{fig:3-6}
        \end{figure}
        \section{联合重建法（SIRT）}
        联合迭代法（SIRT）在ART被提出后不久由Gilbert提出\cite{gilbert1972iterative}，它是ART算法的一
        种并行计算形式。在这种方法中，首先对某个投影$p$所有的像素进行计算之后再对整个重建区域的体素进
        行更新。在更新值加到体素值上之前，需要使用加权值进行加权和归一化。因此SIRT算法收敛速度较慢，
        没有被很广泛地使用。SIRT算法的提出在于使得重建对测量误差不敏感。如上章所述，ART算法每次迭代只
        用到一条射线的射线和（或射线投影）。如果这一射线投影包含误差，则所得的解跟着也引入误差。

        如上文中\ref{eq:2-19}中所示，如果$W$为非奇异矩阵，则其最小二乘法解为
        \begin{equation}
            \bm{x} = (\bm{W}^T\bm{W})^{-1}\bm{W}^T\bm{p}
            \label{eq:3-8}
        \end{equation}
        其中$\bm{W}^T\bm{p}$代表$\bm{p}$的反投影运算。如果把$(\bm{W}^T\bm{W})^{-1}$看作是一个二维滤波器
        ，则式\ref{eq:3-8}就是前面提到的二维滤波反投影。式\ref{eq:3-8}又可以化成
        \begin{equation}
            \bm{x} = \bm{W}^T(\bm{W}\bm{W}^T)^{-1}\bm{p}
            \label{eq:3-9}
        \end{equation}
        式\ref{eq:3-9}中，$(\bm{W}\bm{W}^T)^{-1}$为一维$\rho$滤波器，对$\bm{p}$进行滤波。对式
        \ref{eq:3-9}进行迭代法求解：
        \begin{equation}
            \left\{
                \begin{aligned}
                    & \bm{x}^{0}=\bm{W}^T\bm{p}, \\
                    & \bm{x}^{k+1}=\bm{x}^{k}+\lambda ^{k}(\bm{W}^T\bm{p}-\bm{W}^T\bm{W}\bm{x}^{k})\\
                    & \qquad { } =\bm{x}^{k}+\lambda ^{k}\bm{W}^T(\bm{p}-\bm{W}\bm{x}^{k})\\
                \end{aligned}
                \label{eq:3-10}
                \right.
            \end{equation}
            式\ref{eq:3-10}中，以投影值$\bm{p}$作为反投影的初始值，在第$k+1$次迭代时，利用$k$次的迭代结果
            $\bm{x}^{k}$加上校正值得到$\bm{x}^{k+1}$，校正值与第$k$次估计的误差矢量的反投影
            $\bm{W}^T(\bm{p}-\bm{Wx}^{k})$，因此每个体素的校正值与通过该体素的所有射线的误差值的累加，
            而不只与一条射线相关。因此SIRT的校正过程称为逐点校正（Point by point correction）。这是与ART
            算法的最大的区别，也是SIRT算法能够抑制噪声的根本原因：一些随机误差被通过该体素的所有射线的共
            同贡献平均掉。为了便于迭代计算，SIRT也可改写成：
            \begin{equation}
                x_j^{k+1}=x_j^k + \lambda\sum\limits_{p_i \in P_{\varphi}}\left( p_i -
                \sum\limits_{n=1}^{N}w_{in}x_n^k\right)w_{ij}
                \label{eq:3-11}
            \end{equation}

            图\ref{fig:3-7}是SIRT使用不同的加权模型在$z=100$处的重建结果，重建的参数为：$\lambda=2$（增加
            松弛因子已加快收敛速度），投影角度个数为90，迭代次数20次。
            \begin{figure}[htpb]
                \centering
                \subfigure[使用0,1加权]
                {
                    \label{fig:3-7subfig:1}
                    \includegraphics[width=2in]{3-7-1.jpg} %文件名中不能出现括号
                }
                \subfigure[使用长度加权]
                {
                    \label{fig:3-7subfig:2}
                    \includegraphics[width=2in]{3-7-2.jpg} %文件名中不能出现括号
                }
                \subfigure[线性插值核加权]
                {
                    \label{fig:3-7subfig:3}
                    \includegraphics[width=2in]{3-7-3.jpg} %文件名中不能出现括号
                }
                \caption{SIRT重建结果}
                \label{fig:3-7}
            \end{figure}
            \section{联合代数重建法（SART）}
            1984年，Andersen和Kak 注意到在一次迭代中同一个体素可能被很多光线经过
            \cite{andersen1984simultaneous}。他们证明ART算法中按照光线经过体素的顺序进行体素更新的方法会
            导致重建结果的条状伪影。他们将这些伪影的产生归咎于当时采集图像的噪声。而SIRT算法是在进行一次
            迭代之后对所有的体素一次性更新，这样可以有效地消除条状伪影，但是如上章所示SIRT算法的问题在于
            收敛的速度较慢，需要比较多的迭代次数才能达到较好的效果。Andersen和Kak将ART和SIRT算法的优点进
            行合并就得到SART算法\cite{andersen1984simultaneous}。SART算法不是对投影图像的每个像素（即每条
            光线）进行矫正，而是先计算整个重建区域的投影图像（在角度$\varphi$记为$P_\varphi$），$P_\varphi$中的每个像
            素值对每个体素的矫正值都有贡献，将这些贡献在每个体素上累加就得到了每个体素的更新值。如果简单
            的把矫正项进行相加的话，在投影图像中可能存在的噪声会加入到重建图像中产生伪影，所以在进行更新
            的时候需要进行加权\cite{hj2006}。如图\ref{fig:3-8}为SART投影矩阵计算的示意图。
            \begin{figure}[htpb]
                \centering
                \includegraphics[width=3in]{3-8.jpg} %文件名中不能出现括号
                \caption{SART投影矩阵计算示意图}
                \label{fig:3-8}
            \end{figure}
            SART对重建区域$\bm{x}$的更新方式可以表达为一下形式：
            \begin{equation}
                x_j^{k+1}=x_j^k + \frac{\sum\limits_{p_i \in P_{\varphi}}\left(\lambda \frac{p_i - \sum\limits_{n=1}^{N}w_{in}x_n^k}{\sum\limits_{n=1}^N w_{in}}\right)w_{ij}
            }{\sum\limits_{p_i \in P_{\varphi}}w_{ij}}
            \label{eq:3-12}
        \end{equation}

        式\ref{eq:3-12}与式\ref{eq:3-7}有两个显著不同的地方：1. 对于某个特定的体素$x_j$的矫正项通过计算
        投影图像中相邻像素$p_i$并且用系数$w_{ij}$对像素$p_i$对各体素的影响进行加权。2.尽管以Kaczmarz的方法为
        指导的ART方法需要权值平方之和， SART的思想把ART看作体绘制的逆过程\cite{kaufman1993volume}：体
        数据由一系列累加的体绘制过程（投影过程）以及体绘制的逆过程（反投影）重建
        \cite{mueller1998fast}。之后的研究发现SART能够抑制噪声状的伪影并不是因为基于图像的SART矫正策
        略本身，而是Andersen和Kak用了更好的插值核（例如双线性插值）。但是SART还是有其本身的优势：(1)
        SART更加适合锥束重建；(2)	SART非常适合GPU加速。

        图\ref{fig:3-9}是SART使用不同的加权模型在$z=100$处的重建结果，重建的参数为：$\lambda=0.25$，
        投影角度个数为90，迭代次数5次。
        \begin{figure}[htpb]
            \centering
            \subfigure[使用0,1加权]
            {
                \label{fig:3-9subfig:1}
                \includegraphics[width=2in]{3-9-1.jpg} %文件名中不能出现括号
            }
            \subfigure[使用长度加权]
            {
                \label{fig:3-9subfig:2}
                \includegraphics[width=2in]{3-9-2.jpg} %文件名中不能出现括号
            }
            \subfigure[线性插值核加权]
            {
                \label{fig:3-9subfig:3}
                \includegraphics[width=2in]{3-9-3.jpg} %文件名中不能出现括号
            }
            \caption{SART重建结果}
            \label{fig:3-9}
        \end{figure}
        \section{实验结果分析}
        在前三节中我们已经讨论过经典三种经典迭代重建算法的原理以及实现方法，在本节中我们将通过几个客观评
        价指标来衡量各个算法的优劣\cite{ztg1992ct}。

        在第二章中已经介绍过几种客观评价的方法，本节中我们首先使用Profile曲线法客观评价各个算法的重建图像
        质量，再使用归一化均方根距离测量值$d$（式\ref{eq:2-21}）、归一化平均绝对距离测量值$r$（式
        \ref{eq:2-22}）和图像相似性系数$\varepsilon$（式\ref{eq:2-23}）对各个迭代重建算法的收敛速度和重建质量进行评估。
        \begin{figure}[htpb]
            \centering
            \subfigure[ART]
            {
                \label{fig:3-10subfig:1}
                \includegraphics[width=4.5in]{3-10-1.jpg} %文件名中不能出现括号
            }
            \subfigure[SIRT]
            {
                \label{fig:3-10subfig:2}
                \includegraphics[width=4.5in]{3-10-2.jpg} %文件名中不能出现括号
            }
            \subfigure[SART]
            {
                \label{fig:3-10subfig:3}
                \includegraphics[width=4.5in]{3-10-3.jpg} %文件名中不能出现括号
            }
            \caption{$z=100$切片图像的中心线侧轮廓}
            \label{fig:3-10}
        \end{figure}
        首先为了比较几种重建算法的重建质量，采用Profile曲线法，如图\ref{fig:3-10}所示，选择三种重建算
        法使用线性插值加权的重建结果进行比较，选取$z=100$处的中心线，在直角坐标系中按照对应图像的灰度
        值变化，分别作出灰度变化曲线。从图\ref{fig:3-10}中可以看出，SART算法与其他两种算法相比与标准图
        像拟合度较高，说明SART重建图像的质量较好。

        然后我们再通过$d$、$r$的变化来对三种重建算法的收敛速度进行比较。从表
        \ref{tb:3-1}、\ref{tb:3-2}、\ref{tb:3-3}可以看出，ART和SART重建的收敛速度较快，基本上可以在5次迭
        代的时候得到最优解，而SIRT重建算法的收敛速度较慢，在使用较大的松弛系数的情况下任然需要10次以
        上的迭代才能收敛，也证明了上面的理论研究。同时如图\ref{fig:3-11}所示，随着迭代次数的增加，
        ART重建算法在达到最优结果之后随着迭代次数增加重建结果逐渐偏离。
        SIRT重建算法则收敛速度较慢，在达到最优解之后，重建结果不随着迭代次数增加而显著变差。SART重建
        算法则结合了上述两个算法的有点，不仅能够较快的收敛得到最优解，并且能够在收敛之后重建的结果并
        不随着迭代次数的增加而质量变差。而且从图\ref{fig:3-11subfig:3}中的图像相似性系数$\varepsilon$
        曲线可以看出SART重建算法的结果与标准图像的相似性较高。从图
        \ref{fig:3-11subfig:1}和\ref{fig:3-11subfig:2}中的归一化均方根距离测量值$d$和归一化平均绝对距离
        测量值$r$曲线可以看出SART重建的结果不管是个别点的大误差和整体的小误差都优于其他两种算法的重建
        结果，因此我们在下一章节的文章中主要以SART算法为研究对象进行更加深入的研究。
        \begin{table}[htpb]
            \centering
            \caption{ART重建重建各项指标随迭代次数的变化}
            \begin{tabular}{|c|c|c|c|c|c|}
                \hline
                迭代次数 & 1 & 5 & 10 & 15 & 20 \\
                \hline
                重建误差d& 0.5264 & 0.5006& 0.5727& 0.6305 &0.6859\\
                \hline
                重建误差r& 0.4798 & 0.3830 & 0.4444 & 0.4865&0.5284\\
                \hline
                图像相似性系数$\varepsilon$& 0.8710 &0.9161&0.8878&0.8482&0.8101\\
                \hline
            \end{tabular}
            \label{tb:3-1}
        \end{table}
        \begin{table}[htpb]
            \centering
            \caption{SIRT重建各项指标随迭代次数的变化}
            \begin{tabular}{|c|c|c|c|c|c|}
                \hline
                迭代次数 & 1 & 5 & 10 & 15 & 20 \\
                \hline
                重建误差d& 1.439 & 0.9236 & 0.4238 & 0.4361&0.4455\\
                \hline
                重建误差r& 1.889 & 1.219 & 0.2485 & 0.3758&0.3758\\
                \hline
                图像相似性系数$\varepsilon$& 0.7671 &0.8735&0.9103&0.9181&0.92\\
                \hline
            \end{tabular}
            \label{tb:3-2}
        \end{table}
        \begin{table}[htpb]
            \centering
            \caption{SART重建各项指标随迭代次数的变化}
            \begin{tabular}{|c|c|c|c|c|c|}
                \hline
                迭代次数 & 1 & 5 & 10 & 15 & 20 \\
                \hline
                重建误差d& 0.7176 & 0.3688 & 0.3760 & 0.3852&0.3886\\
                \hline
                重建误差r& 0.7904 & 0.2267 & 0.2586 & 0.2758&0.2796\\
                \hline
                图像相似性系数$\varepsilon$& 0.8616 &0.9345&0.9444&0.9459&0.9448\\
                \hline
            \end{tabular}
            \label{tb:3-3}
        \end{table}
        \begin{figure}[htpb]
            \centering
            \subfigure[归一化均方根距离测量值$d$]
            {
                \label{fig:3-11subfig:1}
                \includegraphics[width=3in]{3-11-1.jpg} %文件名中不能出现括号
            }
            \subfigure[归一化平均绝对距离测量值$r$]
            {
                \label{fig:3-11subfig:2}
                \includegraphics[width=3in]{3-11-2.jpg} %文件名中不能出现括号
            }
            \subfigure[图像相似性系数$\varepsilon$]
            {
                \label{fig:3-11subfig:3}
                \includegraphics[width=3in]{3-11-3.jpg} %文件名中不能出现括号
            }
            \caption{$z=100$切片图像的客观评价分析}
            \label{fig:3-11}
        \end{figure}
        \section{本章小结}
        在本章中，主要以二维情况为例讨论了三种经典迭代算法的各自特点，最后以4个客观评价指标：Profiler曲线
        、归一化均方根距离测量值$d$、归一化平均绝对距离测量值$r$、图像相似性系数$\varepsilon$对重建结果做
        评价，进一步比较各个重建算法的优劣进行比较。

        本章首先讨论了投影矩阵的计算方法，介绍了4种最常用的投影矩阵计算方法以及各自特点，其中着重介绍了基
        于线性插值的投影矩阵计算，因为在之后的迭代算法的研究中都是以此方法计算得到的投影矩阵进行试验的。

        然后分别从理论以及从实现的角度介绍了ART、SIRT、SART三种重建算法，并且以0、1加权、长度加权、
        和线性插值加权计算的投影矩阵用于三种重建算法进行重建试验。

        最后一节中以4个客观评价指标对三种重建算法的重建结果进行比较，并且比较了随着迭代次数的增加，各个客
        观指标的变化情况，最后发现SART不管是从迭代收敛速度还是重建的质量都是三者中最佳的，为下一步以SART
        算法为主要研究对象做铺垫。
        \chapter{SART算法的研究}
        上章中讨论了二维平行射线情况下的三种经典迭代重建算法各自的特点，最终得到SART是一种较为优秀的重建
        算法的结论。本章主要讨论的内容是将SART重建算法推广到三维锥束情况下的各项重建指标。

        在锥束重建的情况下，SART算法相对于ART具有另一项优点：由于ART算法是逐线更新的重建过程，在平行的情
        况下，各个插值核所覆盖的重建区域是相同的。但是在锥束情况下，插值核在距离射线源近的地方覆盖的相对于
        距离射线源较远的地方更加密集，距离射线源更近的的采样点的加权值会比较大，从而会导致反投影时能量会
        更加集中在距离射线源更近的地方，如图\ref{fig:4-1}所示，最后影响重建图像的质量。要使得锥束ART重建
        算法避免这个问题，需要使用可变的线性插值模型，插值核的大小随着射线上的采样点距离射线源的远近改变
        ，这样不可避免地导致计算复杂度增大。
        \begin{figure}[htpb]
            \centering
            \subfigure[锥束反投影侧视图]
            {
                \label{fig:4-1subfig:1}
                \includegraphics[width=2.5in]{4-1-1.jpg} %文件名中不能出现括号
            }
            \hspace{2cm}
            \subfigure[平行反投影侧视图]
            {
                \label{fig:4-1subfig:2}
                \includegraphics[width=2.5in]{4-1-2.jpg} %文件名中不能出现括号
            }
            \caption{锥束ART重建反投影侧视图}
            \label{fig:4-1}
        \end{figure}
        而在SART重建算法进行锥束重建时，它所使用的是逐点更新的策略，即是不管采样点距离射线源远近，最后投
        影反投影加权时都是跟在同一个体素中的采样点加权，同一体素内部的采样点距离射线源基本相同，所以加权
        值也在相同的尺度上，最后加权之后不会导致能量分布不均（如图\ref{fig:4-2}）所示，因此SART只需要固定
        大小的插值核即可进行重建。具体证明过程可参考文献\cite{mueller1998fast}。
        \begin{figure}[htpb]
            \centering
            \includegraphics[width=2.5in]{4-2.jpg} %文件名中不能出现括号
            \caption{锥束SART重建反投影侧视图}
            \label{fig:4-2}
        \end{figure}
        基于上一章节的结论以及本章之前的讨论，在本章余下的内容中都以SART重建算法进行讨论。
        \section{基于体素驱动的SART}
        在展开讨论之前首先介绍几个术语。在文献中\cite{crawfis1993texture}中，Crawfis和Max提出我们可以把插
        值核插值核想象成一个具有纹理的多边形放置到重建空间中。这里的纹理通过插值核函数计算的。在下文中我
        们将置于重建空间中的插值核多边形简称为插值核，而将插值核多边形在成像平面上的投影称为插值核图像。
        由于当一条射线射入重建区域时，对投影图像上的每个像素值的累加值都相等，所以当把插值核多边形放置到
        重建区域的时候需要放置于与对应射线方向正交的位置。线性积分就通过寻找射线与插值核多边形中相交的点
        的在插值核多边形中索引值获得的。为了加速计算的过程，我们通常根据确定的插值核预先计算一个积分查找
        表，这个查找表的索引值与具体插值计算方法相关，根据索引值查找得到当前计算的射线与体素的加权系数，
        如图\ref{fig:4-3}所示，这是一个二维透视投影积分查找表的积分值分布图。当射线穿越插值核时可以直接通
        过查找表来查找积分值，因此为了splatting精确，需要把二维的插值核映射到射线垂直穿过插值核时的像素值
        。只有这样真正的光线插值核积分与预计算的查找积分表相匹配。但是在实际运用中为了加快投影矩阵计算的
        速度，可以将上述的条件适当放宽，做一定的近似。Westover在文献\cite{westover1989interactive}中提出
        的体素驱动splatting与上文所述的有三处不同：

        （1）在放置插值核多边形计算它在投影平面上的投影范围时并不是将多边形放置到与中心射线垂直的位置，而是
        将它放置得与投影平面平行；

        （2）当计算插值核到投影平面像素的映射时以平行投影取代透视投影（如图\ref{fig:4-4subfig:1}所示）；

        （3）插值核多边形不因为需要与穿越插值核多边形的射线垂直而旋转（如图\ref{fig:4-4subfig:2}所示）。
        \begin{figure}[htpb]
            \centering
            \includegraphics[width=4in]{4-3.jpg} %文件名中不能出现括号
            \caption{插值核查找表积分分布图}
            \label{fig:4-3}
        \end{figure}
        \begin{figure}[htpb]
            \centering
            \subfigure[平行投影取代透视投影]
            {
                \label{fig:4-4subfig:1}
                \includegraphics[width=4in]{4-4-1.jpg} %文件名中不能出现括号
            }
            \hspace{2cm}
            \subfigure[不旋转插值核多边形]
            {
                \label{fig:4-4subfig:2}
                \includegraphics[width=4in]{4-4-2.jpg} %文件名中不能出现括号
            }
            \caption{Westover两种简化情况}
            \label{fig:4-4}
        \end{figure}

        其中第三点的影响是非常细微的，之前的两点相对比较更加重要。第一点近似计算插值核在投影平面上
        的投影大小时比实际的情况略小。比如锥角为$30^o$重建区域大小为$128^3$的情况下，插值核长度的精确值与
        近似值之间的相对误差率最大时1.15，最大的绝对误差是1.46个像素。最大的误差发生在距离发射源较近在锥
        束边界附近的那些体素上。这个误差导致了这些体素在投影平面上覆盖的区域比实际情况略小。第二种近似具
        有类似的影响。在这种情况下插值核积分查找表到投影平面像素的映射会略微被压缩。同样的离发射源较近且
        在锥束边界的体素受影响较大。

        如图\ref{fig:4-5}所示，这是对透视体素驱动splatting的精确解法。为了演示简便，图中只显示了2D的情
        况。我们首先把坐标系的原点定位在射线源的位置。在重建的过程中，对于不同角度的投影我们首先将重建区
        域旋转与投影对应的角度。将插值核多边形放置于体素$v_{x,y,z}$处，并且与射线源到体素$v_{x,y,z}$中心的
        向量正交。注意到这样的策略只有在经过中央的直线的线性积分是精确的，其他的射线经过体素的核函数时的
        方向会与放置在重建空间中的二维（图\ref{fig:4-5}中是一维）插值核多边形给的方向会有略微不同。因此
        Westover提出的第三个近似的误差还存在，但是这个误差即是对于离射线源近的体素也不到$1\%$，因此可以
        忽略。
        \begin{figure}[htpb]
            \centering
            \includegraphics[width=4.5in]{4-5.jpg} %文件名中不能出现括号
            \caption{体素驱动投影矩阵计算示意图}
            \label{fig:4-5}
        \end{figure}

        插值核多边形所在的平面等式的参数由归一化的中央射线确定（即射线源到体素$v_{x,y,z}$中心的向量）。从这
        个等式中我们可以得到在同一平面的两个正交向量：$u$和$w$（图\ref{fig:4-3}中只显示了$u$）。这里使用
        $u$和$w$是为了将他们投影到投影平面的两个轴上。通过$u$和$w$可以在投影平面上找到插值核多边形的四个顶
        点在空间坐标系中的坐标（图\ref{fig:4-3}显示其中的两个：$V_{Right}(v_{x,y})$和$V_{Left}(v_{x,y})$
        ）。这四个顶点被透射投影到投影平面上，在投影平面上产生里一个插值核图像的矩形范围（图\ref{fig:4-3}显
        示其中的两个：$Ext_{Right}(v_{x,y})$和$Ext_{Left}(v_{x,y})$）。将经过插值核多边形的射线与多边形的交
        点以参数的形式表达，则可制定一套递增的策略将插值核图像中的像素值与插值核积分查找表的索引值联系起来，这样
        可以减少一定的计算复杂度。

        将插值核多边形映射到投影平面和建立递增的将投影图像的像素映射到插值核查找表的策略都具有相当的计算复杂
        度：大约有100个乘法、加法和除法，以及两个开方运算。如果不使用递增的策略可以用来加速邻接的体素的映
        射，这样的映射操作通常时间复杂度为$O(N)=O(n^3)$。在本文的实现中，透视投影的运算时间复杂度大约
        是平行投影的两倍。
        \section{基于射线驱动的SART}
        在上一节中透视体素驱动splatting可以产生精确重建的投影矩阵，但是基于这种方法的锥束重建非常耗时。
        而基于射线驱动的方法对于透视变换的非线性性不敏感，因此基于射线驱动可能可以更好地利用锥束透视投影
        的特性。Hanson和Wecksung在文献\cite{hanson1985local}中将二维的射线驱动算法扩展到三维空间中去。这
        个算法对使用固定插值核的锥束SART算法和使用可变插值核的锥束ART算法都适用。

        在射线驱动算法中，体素不再同时对所有的射线在投影平面的像素进行累加，而是每个像素对于的射线分别累
        加，这样相比体素驱动更加适合ART算法的实现。以下就是射线驱动实现的过程：首先把重建区域划分成相互平
        行且与投影平面平行的二维平面。当一条射线射入重建区域时，它与每个二维平面相交并且决定在这个平
        面内被射线经过的体素。如图\ref{fig:4-6subfig:1}所示，与投影平面像素$p_i$对于的射线射入重建区域，与平面
        $x_s$相交于$y=y(i,x_s)$。在平面$x_s$中的与射线相交的体素核由区间这里$\alpha$是射线的斜率。
        \begin{equation}
            y_{Right}(i,x_s)=y(i,x_s)+\frac{extent_{kernel}}{cos(\alpha)}
            \label{eq:4-1}
        \end{equation}
        而$y_{Left}(i,x_s)$类似。在确定有效的体素的区间$[Ceil(y_{Left}(i,x_s)),Floor(y_{Right}(i,x_s)]$之
        后在平面$x_s$上且在区间范围内的体素被称为有效体素，接下去需要计算有效体素在插值核查询表中的索引值
        。如图\ref{fig:4-6subfig:2}所示，对于一个位于坐标$(x_v,y_v)$处的体素$v$，它在插值核查找表中的索引
        值可以由穿越体素中心的且射线平行的直线与射线之间的距离$dr$确定。$y_{Right}(i,x_s)$由公式
        \ref{eq:4-1}确定。$dr$可由式\ref{eq:4-2}确定。
        \begin{equation}
            dr=a\cdot x_s+b\cdot y(i,x_s)-a\cdot x_s - b\cdot y_v = (b\cdot (y(i,x_s) -y_v)
            \label{eq:4-2}
        \end{equation}
        这里$a$和$b$是射线的直线方程$a\cdot x_s+b\cdot y(i,x_s) = 0$的参数，也可由射线的归一化向量求得。
        在射线在层之间遍历时可以在前一层的$y_{Left}(i,x)$，$y_{Right}(i,x)$和$dr$的基础上进行累加，无需
        重新进行计算。
        \begin{figure}[htpb]
            \centering
            \subfigure[确定平面上的有效体素]
            {
                \label{fig:4-6subfig:1}
                \includegraphics[width=4in]{4-6-1.jpg} %文件名中不能出现括号
            }
            \hspace{2cm}
            \subfigure[计算有效体素的索引值]
            {
                \label{fig:4-6subfig:2}
                \includegraphics[width=4in]{4-6-2.jpg} %文件名中不能出现括号
            }
            \caption{光线驱动计算投影矩阵示意图}
            \label{fig:4-6}
        \end{figure}
        在三维的情况中，一条直线需要两个平面确定。在这里，一条三维的射线由两个正交的平面确定。其中一个平
        面的法向量由直线的向量与和投影平面上的任一坐标轴叉乘所得。另一平面的法向量由射线向量与第一个二维
        切片平面的法向量叉乘所得。这两个平面相互垂直而在这两个平面上的插值核多边形也相互垂直。水平的平面
        与插值核多边形的平面相交再使用公式\ref{eq:4-2}计算的水平方向的索引值$dr_{row}$。以垂直方向的平面
        与插值核多边形相交可计算得垂直方向的索引值$dr_{col}$，根据这两个索引值在索引查找表中的值即是该射
        线对于体素$v$的线性积分值。

        实现上述的过程，需要三层循环。最外层循环初始化射线射入到重建区域中。中间层循环控制射线在分层之后
        的重建区域内遍历。最里层循环再确定的范围之内的体素对与当前射线的贡献。对于透视投影，两个平面的等
        式的参数每条射线都需要计算一遍。一次需要多50次加法、乘法和除法，以及三次开方运算。一条射线在层与
        层之间遍历决定有效体素以及计算插值索引表中的值的运算量与在体素驱动中旋转插值核再splatting到投影屏
        幕上确定插值索引表中的索引值的计算量相当。射线驱动将splatting算法按照体素遍历改成按照光线遍历，因
        此最外层算法的时间复杂度为$O(n^2)$。这样在透视投影的情况下，为透视投影调整而需要额外的计算比基于
        体素驱动的算法少了一个数量级（体素驱动的复杂度为$O(n^3)$）。而且更重要的一点是光线驱动中不需要影响
        精度的估计。事实上，此算法从设计上就避免了体素驱动方法中的索引误差。
        \section{锥束SART重建结果分析}
        在本章中，我们首先将对两种投影矩阵的计算方法进行比较。在实现基于体素驱动的投影计算时使用了
        Westover提出的三种近似情况以提高算法效率，而射线驱动的算法直接按照上节中所提出的算法进行实现，首
        先先比较两种投影矩阵方法的计算效率根据表\ref{tb:4-1}、\ref{tb:4-2}展示的实验结果,从总体上看基于射
        线驱动的投影矩阵计算效率高出不少。而且当重建区域体素越多时，这种优势越明显。所以基于射线驱动的
        SART算法更加适合后面的较高分辨率的重
        建场景。因此在以下的实验中，以射线驱动计算的投影矩阵进行SART重建实验。
        \begin{table}[htpb]
            \centering
            \caption{基于体素驱动在各个分辨率下计算投影矩阵所需的时间（单位：ms）}
            \begin{tabular}{|l|*{4}{c|}}\hline
                \backslashbox{重建区域}{投影图像} & 64$\times$64 & 128$\times$128 & 256$\times$256 \\\hline
                64$\times$64$\times$ 64 &40.3252&120.532&531.374 \\\hline
                128$\times$128$\times$128&314.535&1018.14&4058.92\\\hline
                256$\times$256$\times$256&2547.73&8247.01&32877.2\\\hline
            \end{tabular}
            \label{tb:4-1}
        \end{table}
        \begin{table}[htpb]
            \centering
            \caption{基于射线驱动在各个分辨率下计算投影矩阵所需的时间（单位：ms）}
            \begin{tabular}{|l|*{4}{c|}}\hline
                \backslashbox{重建区域}{投影图像} & 64$\times$64 & 128$\times$128 & 256$\times$256 \\\hline
                64$\times$64$\times$ 64 &73.1712& 140.338 & 453.694 \\\hline
                128$\times$128$\times$128&139.284& 554.86 & 2012.22\\\hline
                256$\times$256$\times$256&266.891 &1085.36 &4329.28 \\\hline
            \end{tabular}
            \label{tb:4-2}
        \end{table}

        在投影矩阵计算完成之后即可按照各个投影角度的投影矩阵进行重建。锥束SART重建一次迭代的流程图如
        \ref{fig:4-7}所示。当完成一次迭代之后，以一次迭代的结果为初始值进行下一次迭代，直到迭代次数达到预
        设值。SART再各个投影个数下的重建结果如图\ref{fig:4-8}所示。
        \begin{figure}[htpb]
            \centering
            \includegraphics[width=2in]{4-7.jpg} %文件名中不能出现括号
            \caption{迭代重建一次迭代流程图}
            \label{fig:4-7}
        \end{figure}
        \begin{figure}[htpb]
            \subfigure[]
            {
                \label{fig:4-8subfig:1}
                \begin{minipage}{0.3\textwidth}
                    \includegraphics[width=2in]{4-8-1.jpg}\\
                    \vspace{0.1in}
                    \includegraphics[width=2in]{4-8-4.png}\\
                    \vspace{0.1in}
                    \includegraphics[width=2in]{4-8-7.png}
                \end{minipage}%
            }
            \hspace{0.1in}
            \subfigure[]
            {
                \label{fig:4-8subfig:2}
                \begin{minipage}{0.3\textwidth}
                    \includegraphics[width=2in]{4-8-2.jpg} \\
                    \vspace{0.1in}
                    \includegraphics[width=2in]{4-8-5.png} \\
                    \vspace{0.1in}
                    \includegraphics[width=2in]{4-8-8.png}
                \end{minipage}
            }
            \hspace{0.1in}
            \subfigure[]
            {
                \label{fig:4-8subfig:3}
                \begin{minipage}{0.3\textwidth}
                    \includegraphics[width=2in]{4-8-3.jpg} \\
                    \vspace{0.1in}
                    \includegraphics[width=2in]{4-8-6.png} \\
                    \vspace{0.1in}
                    \includegraphics[width=2in]{4-8-9.png}
                \end{minipage}
            }
            \caption{分别30、60,、90个角度时(a)z=100(b)y=134(c)x=134迭代1次的重建结果}
            \label{fig:4-8}
        \end{figure}
        在各个角度重建所需时间如表\ref{tb:4-3}所示（实验平台为Intel(R) Core(TM) Q8400(2.66GHz)，处理器4GB
        内存，投影分辨率为256$\times$256，重建区域分辨率为$256\times 256 \times 256$），由表可见SART重建
        的速度非常慢，而且随着重建区域分辨率以及投影射线个数的增大，重建的时间呈指数级上涨而重建所需的内
        存空间也会大大的增加。在实际运用中，当投影设备分辨率较高或者要求重建区域需要较高分辨率时，SART等
        迭代重建算法将无法满足实际工程应用的需要。
        \begin{table}[htpb]
            \centering
            \caption{基于射线驱动在各投影个数下重建所需的时间（单位：s）}
            \begin{tabular}{|l|*{2}{c|}}\hline
                投影个数 & 重建时间\\\hline
                30 &233.562\\\hline
                60 &452.374\\\hline
                90 &686.251\\\hline
            \end{tabular}
            \label{tb:4-3}
        \end{table}
        \section{本章小结}
        本章主要介绍了锥束SART重建的实现过程，主要介绍在锥束情况下两种投影矩阵的计算方法：体素驱动和射线
        驱动。基于体素驱动的投影矩阵计算方法中，利用Westover提出的三种近似方法以降低计算的时间复杂度以及
        实现算法的逻辑。与此相比基于射线驱动的投影矩阵计算方法是完全精确的计算方法，不需要进行近似精简，
        而且在处理透视投影时更加简单。随后分别对两种投影矩阵计算方法进行实验，实验结果表明基于射线驱动的
        投影矩阵计算方法更加快速精确。之后使用基于射线驱动的投影矩阵进行SART重建的实验，实验结果表明SART
        重建的结果具有较好的鲁棒性，且能够利用迭代算法在较少投影的情况下获得比较好的重建结果。但是与解析
        类型的重建方法相比，迭代重建的速度还是相当慢而且在重建区域分辨率较高时对内存的要求很高，所以在实
        时性要求较高的应用中无法达到应用要求。
        \chapter{SART算法的CUDA加速}
        在上一章节中介绍了SART重建的基本原理和实现，但是由于迭代算法本身的特点，使用SART进行重建的速度非
        常慢，而且对内存具有较高的要求，因此并不适合实际工程的应用，在本章中将利用CUDA开发环境对进行一定
        优化过之后的SART重建算法进行GPU加速，使得利用SART算法进行迭代重建的具有实际运用的可行性。在本章中
        首先介绍了CUDA开发环境的基本概念，之后再对SART算法进行一定的改进以适应CUDA程序的开发要求，最后将
        串行的CPU代码改成并行的CUDA代码，并且按照CUDA的一些特性，对并行计算的任务划分进行一定的优化，以达到尽
        可能大的加速比。
        \section{基于CUDA的GPU通用计算技术}
        \subsection{GPU通用计算的发展}
        GPU（Graphic Processing Unit）全程图形处理器，在上世纪八九十年代，GPU主要的主要用作工作站的计算机
        图形特别是三维图形的渲染，此时图形处理流水线的大部分功能尚由CPU处理，随着GPU功能的不断完善CPU的
        流水线功能不断得向GPU转移，同时加入了各种高级功能，并且形成了一套标准的图形学开发库——
        OpenGL\cite{unit2004}。到了
        90年代末期，随着消费电子（如电子游戏）需求的不断旺盛，GPU的迎来了发展的黄金时期，从98年第一款桌面
        独立图像处理3DFX Voodoo3发布直到目前为止，平均一年就有新一代图形卡诞生。随着GPU渲染能力的不断提升
        ，利用GPU进行大规模并行计算的潜力渐渐受到科研人员的关注。

        GPU通用计算（GPGPU，General-Purpose Computing on Graphics Processing Units）主要采用CPU+GPU异构模
        式，CPU主要处理事务、逻辑事件响应、为GPU准备数据等串行计算，GPU进行计算密集型的大规模并行计算。与
        CPU相比，GPU具有相当数量的计算核心以及更大的数据带宽在进行大规模数据计算时具有较大的优势。传统的
        GPU构架如图\ref{fig:5-1}所示，GPU是由一系列流水管线组成的。其中顶点投影以及纹理装配过程是可编程的
        ，GPU通用计算最原始的实现就是利用这两个步骤的可编程性，将数据打包成顶点数据或者纹理数据再通过流水
        管线进行并行计算，所以一开始的GPU通用计算都是基于汇编或者图形学API的，最典型的基于图形学API是基于
        OpenGL扩展的GLSL（Graphic Library Shader Language）\cite{rost2005opengl}。
        \begin{figure}[htpb]
            \centering
            \includegraphics[width=4in]{5-1.jpg} %文件名中不能出现括号
            \caption{图形学渲染管线示意图}
            \label{fig:5-1}
        \end{figure}

        CUDA是由NVIDIA公司于2007年6月发布的一套体系。这套体系包括软件和硬件两部分，主要用于将GPU作为数据
        并行计算设备。到目前为止，CUDA的版本已经进行了4次大的版本升级，功能不断完善\cite{zs2009gpu}。CUDA相对于传统的构架
        有三个显著的进步：一是采用统一处理构架（Unified Graphics and Computing Architecture）可将一个图形
        处理单元看做一个类似CPU的计算核心，二是引入片内共享存储器，支持随机读写和线程通信，三是提供一套
        基于ANIS C标准扩展的编程接口，使得普通C++程序员更快掌握GPU通用计算编程。
        \subsection{CUDA编程模型}
        一个CUDA程序包括主机代码与设备代码，NVIDIA C 编译器（nvcc）在编译阶段把两部分代码分开。主机端代码
        是普通的ANIS C代码，标准的C编译器即可编译成普通的主机端线程。而用ANIS C扩展标记写成的设备端代码需
        要特殊的处理，这些代码主要定义数据并行计算的函数以及与之相关的数据结构，这些函数被称为$kernel$。
        设备代码将进一步被nvcc编译成可在GPU上执行的代码。CUDA程序运行时包括若干个不同的运行状态：在主机端
        执行（如CPU）和在设备端执行（如GPU）\cite{hwu2009programming}。在一个系统中可以包含一个主机与若干
        个设备。在最新发布的CUDA标准中提供了统一虚拟寻址（UVA，Unified Virtual Addressing）的特性，可以将
        主机端存储器空间与设备端存储器空间整合
        为统一的寻址空间。图\ref{fig:5-2}为CPU+GPU异构编程模型示意图，在主机端执行的状态时程序主要串行执
        行主要完成在$kernel$启动之前数据准备与设备初始化操作，在设备端状态时程序执行大规模的并行计算。
        \begin{figure}[htpb]
            \centering
            \includegraphics[width=4in]{5-2.jpg} %文件名中不能出现括号
            \caption{CPU+GPU异构编程模型}
            \label{fig:5-2}
        \end{figure}

        $kernel$函数启动时将产生大量的并行执行的线程，每个$kernel$对应于一个$Grid$，每个$Grid$由若干个线
        程块$block$组成，而若干个$thread$组成一个$block$。$block$在$Grid$中的维数与
        $thread$在$block$中的维数在$kernel$启动时确定，并且在整个$kernel$运行期间不变。这种线程的组织方式
        与支持CUDA的图形卡的硬件结构有关，支持CUDA的GPU主要由若干个SP（Stream Processor）组成，每个SP对应
        执行一个线程，但是SP本身不带任务调度的硬件，它依赖于上层硬件SM（Stream Multiprocessor）每个SM包含
        有8个SP以及进行任务调度的硬件，这种结构被称为SIMT（Single Instruction，Multiple Thread，单指令多
        线程）模型。在硬件执行时，每个$block$会绑定到一个SM上，并且被划分成若干个
        $warp$，每个$warp$包含最多32个$thread$，当一个$warp$执行同步或者显存访问等延时较高的操作时，SM可以切
        换到另外一个$warp$执行，这样可有效地隐藏访问延时，提高GPU的执行效率。CUDA的这种结构设计具有较好的
        可伸缩性，一个在具有较少SP的GPU上设计的CUDA程序换到更好的GPU上执行时，无需修改代码即可得到相应
        的性能提升。
        \subsection{CUDA存储模型}
        在CUDA中，设备的内存空间主要是指显卡上所载的显存空间，一般在1GB到4GB之间。主机端代码能够通过CUDA
        的API对显存直接进行读写。但是CUDA编程时还需要其他存储空间的配合才能发挥最大的效能，图
        \ref{fig:5-3}表示的是CUDA的存储模型（双箭头表示访问权限关系）。每个$thread$都有自己的寄存器（
        Register）和局部存储器（Local Memory），但是每个$thread$所占的寄存器个数不定，因为SM上寄存器的个
        数是固定的，寄存器将平均分配到与SM绑定的$block$中的$thread$中，如果每个$thread$所占寄存器过多时
        $warp$同时运行的$thread$个数将减小，所以如何合理利用寄存器资源是CUDA程序优化的一个重要方面。每个
        $block$都有一块共享存储器（Shared Memory），共享存储器可以被同一个$block$中的所有$thread$访问，因
        此共享存储器是$thread$之间进行通信的重要手段，而且共享存储器的访存速度并全局存储器的访存速度快几
        个数量级，与寄存器访存速度相当。所有$thread$都能够访问全局存储器（Global Memory），在$thread$访问
        全局存储器时CUDA可以自动进行访存的归约，如图\ref{fig:5-4}所示，当相邻的$thread$访问连续的全局存储
        器区域时，本来各个$thread$各自的访问内存的指令会归约成一条访存指令（具体归约的能力与GPU的计算能力
        相关），并且连续读取一整块显存区域的速度会比各个$thread$分别随机读取速度快很多。CUDA访存归约的特
        性配合$block$中的共享存储器可以对很多应用情况进行访存优化：对并行任务划分时尽量以显存中的布局为依
        据进行$block,thread$的划分，在$thread$的初始阶段将所需的数据从全局存储器中读取到共享存储器之中，
        之后使用数据时直接从共享存储器中读取，程序速度将大大提高。除了上诉三种存储器类型之外，还有两种可
        被所有$thread$访问的只读存储器：常数存储器（Constant Memory）和纹理存储器（Texture Memory）。常数存
        储器的容量较小，常被用于存储经常会访问的数据结构，而纹理存储器主要加入了纹理缓存再把一部分全局存
        储器划分为纹理存储器的存储区域，在最新的CUDA构架中，对于全局存储器也加入了访问缓存的支持，所以纹
        理存储器访问速度的优势不在，相对于全局存储器只剩下在利用索引访存时自动进行越界检测和自动进行插值
        的功能优势。
        \begin{figure}[htpb]
            \centering
            \includegraphics[width=4in]{5-3.jpg} %文件名中不能出现括号
            \caption{CUDA存储模型}
            \label{fig:5-3}
        \end{figure}
        \begin{figure}[htpb]
            \centering
            \includegraphics[width=3in]{5-4.jpg} %文件名中不能出现括号
            \caption{CUDA访存归约示意图}
            \label{fig:5-4}
        \end{figure}
        \section{射线驱动的SART改进算法}
        上一小节中介绍了利用CUDA进行GPU并行加速的理论基础和方法，但是在将SART算法进行CUDA加速时需要进行一
        定的改进以使之更加适于改成并行化算法，本节中将主要介绍基于射线驱动的SART算法的改进算法。

        在对SART算法进行并行化加速时主要遇到两个问题：一是在每次旋转到新的角度投影时计算投影矩阵的算法中
        含有逻辑判断较多，不适合并行计算，二是当重建分辨率较高时，对存储容量要求较高，CPU版本运行的代码主
        要使用系统内存，成本较低。当需要在图形卡上运行时，显存的成本较高。传统迭代重建方法的CUDA加速方案
        中，将投影矩阵的计算放在CPU端进行，在重建时将投影矩阵传到显存后在GPU上进行并行计算。此方案在投影
        图像和重建数据分辨率较高时会引起以下问题：

        1） 投影矩阵的计算非常耗时，从而影响整体重建加速的性能。

        2） 投影矩阵较大，导致显存不足，无法进行较高分辨率的重建。

        为了解决以上问题，本文提出了一种新的重建方案。对于任何投影系统的旋转都是刚体变换，所以可以将这种
        旋转映射成被重建区域的旋转。基于这个原理对这个重建过程进行了一定的改进如图\ref{fig:5-5}所示。这样
        的改进具有以下的优势：
        \begin{figure}[htpb]
            \centering
            \includegraphics[width=3in]{5-5.jpg} %文件名中不能出现括号
            \caption{改进后SART重建算法流程}
            \label{fig:5-5}
        \end{figure}

        1）	避免了对每个角度都计算投影矩阵，计算投影矩阵的过程不适合改成并行算法加速，而旋转重建数据
        的过程非常适合并行加速，并且CT旋转可看做刚体变换，旋转矩阵的计算速度相比投影矩阵的计算速度具有明
        显优势。

        2）	当初始角度为0时，由于投影矩阵的对称性，可只保存1/4大小的投影矩阵，其他都可由对称性获得。图
        \ref{fig:5-6}展示射线经过体素某个平面时计算投影矩阵参数时的情况，由于对称性投影矩阵只需要计算和保
        存$\omega_{i,j}$的内容，当投影和反投影需要访问
        $\omega_{i^{'},j},\omega_{i,j^{'}},\omega_{i^{'},j^{'}}$时可以通过对称性访问$\omega_{i,j}$的值即
        可，所以投影矩阵只需要保存左下角1/4大小，其余任何一个条光线对于的矩阵行都可通过对称性计算所得，这样
        不但可以减少投影矩阵的计算时间，也能节约内存空间。这点在高分辨率重建的过程中非常得重要。
        \begin{figure}[htpb]
            \centering
            \includegraphics[width=3in]{5-6.jpg} %文件名中不能出现括号
            \caption{利用对称性分块计算投影矩阵示意图}
            \label{fig:5-6}
        \end{figure}

        对于上述第二点，可以对SART重建的内存使用情况进行分析以证明此举的意义。在SART重建的过程中，主要的显存占用情况包括几个部分：

        1）	投影数据；

        2）	重建数据；

        3）	投影矩阵；

        4）	SART重建过程中的临时数据。

        我们以180个角度$256 \times 256$的投影数据，$256\times 256\times 256$的重建数据为例，进行分析。

	    1）投影数据：假设以16位数据保存，则$M_p=256 \times 256 \times 180 \times 2= 23592960 B=22.5MB$；

        2）重建数据：假设以32位float类型保存，则$M_r=256 \times 256 \times 256 \times4= 67108864 B=64MB$；

        3）投影矩阵：假设用int型保存体素编号（$256\times 256\times 256$投影的话24位数据即可保存，考虑到扩展性这里用36位的结构
        体保存，每一维用12位表示），对应的参数ω由一个float表示，每条光线最多经过的体素个数为307（经实验确定为
        体素宽度乘以1.2），则$M_w=256 \times 256\times 307 \times 12= 241434624B =230MB$；

        4）SART重建过程中的临时数据：主要是在更新过程中的临时数据，包括每个体素被穿过的光线编号以及对应的
        系数值$\omega_{ij}$（假设每个体素最多被5条光线穿过） ，以及实际经过每个体素的光线条数；$M_t=256
        \times 256 \times 256 \times (5\times 4+8)= 469762048 B=44MB$；

        所以SART重建最基本的内存消耗为$M=M_p+M_r+M_w+M_t= 764M$在实验过程中，整个程序的内存的使用超过了1
        G，在如今的PC上这个内存占有量并不大，但是在显卡上已经无法满足需要了。所以需要利用对称性对数据进行
        压缩，压缩后$M_w=57M$，$M_t=164M$，从而让在GPU上进行较高分辨率的重建成为了可能。

        表\ref{tb:5-1}表示原始SART重建算法与改进后的SART重建算法各项评价指标在各个投影个数重建的值，由表
        可看出改进后算法的重建结果与原始算法差别极小，主要误差由于选择时插值计算引入，而压缩投影矩阵的情
        况理论上跟不压缩重建的结果是完全一致的。
        \begin{table}[htpb]
            \centering
            \caption{改进前后SART重建结果各项指标的比较}
            \begin{tabular}{|l|*{4}{c|}}\hline
                \backslashbox{评价指标}{投影个数} & 30 & 60 & 90 \\
                \hline
                重建误差d& 0.0691 & 0.0792 & 0.0760 \\
                \hline
                重建误差r& 0.0442 & 0.0482 & 0.0453 \\
                \hline
                图像相似性系数$\varepsilon$& 0.9980 &0.9972&0.9973\\
                \hline
            \end{tabular}
            \label{tb:5-1}
        \end{table}
        表\ref{tb:5-2}展示了SART重建算法改进前以及改进后压缩与不压缩的投影矩阵情况下各投影个数时的重建时
        间的比较，可看出改进后的SART算法在速度上比改进前具有较大优势，这是因为矩阵旋转的速度比计算一个特
        定角度的投影矩阵速度快很多，而在压缩投影矩阵后由于只需要计算1/4数量的射线所对应的投影矩阵，虽然在
        访问数据时需要一系列的镜像操作，但实际速度比不压缩投影矩阵时仍具有一定速度上优势。
        \begin{table}[htpb]
            \centering
            \caption{改进前后SART重建速度比较（单位：s）}
            \begin{tabular}{|l|*{4}{c|}}\hline
                \backslashbox{重建所需时间}{投影个数} & 30 & 60 & 90 \\
                \hline
                原始SART&233.562&452.374 & 686.251 \\
                \hline
                改进SART不压缩投影矩阵& 73.7152&142.033&218.419\\
                \hline
                改进SART压缩投影矩阵& 69.7640&140.497&215.954\\
                \hline
            \end{tabular}
            \label{tb:5-2}
        \end{table}

        通过以上的实验表明经过改进的SART算法在存储空间和计算速度上较原始SART重建算法都有一定的优势，但是
        这个改进最大的优势在于它更容易被改成并行化程序，首先旋转矩阵的操作是一个非常适合并行化的操作，而
        计算投影矩阵的操作由于矩阵维数的不确定性以及判断逻辑负责难移并行化。所以在下一节中将介绍利用CUDA
        开发环境讲改进后的SART算法改造成在GPU上运行的并行程序。
        \section{SART算法CUDA加速}
        \subsection{SART的CUDA加速算法的任务划分}
        在将SART重建算法改成并行算法的过程中，主要包含对四个部分算法的改进：旋转重建数据，正投影过程，反
        投影过程，逆旋转重建数据，下面将分别对各个部分并行化过程进行阐述。

        旋转与你旋转重建数据的过程是相似的过程，所以讲这两个算法一并进行介绍。在旋转之前首先需要求一个旋
        转矩阵，因为旋转过程使用的是双线性插值，所以矩阵的内容包括被旋转体素上下左右4个体素的值（如图
        \ref{fig:5-7}所示），由于每个被旋转体素的旋转后的值是相互独立的，所以非常适合并行计算，只需使用一
        个$thread$计算相对于的旋转体素的计算即可。在旋转时需要经常访问邻近体素的信息，所以可以在每个线程
        开始时将邻近体素读取到共享内存中（Shared Memory），图\ref{fig:5-8}展示的就是一个$block$所对应的体
        素，将在每个$block$开始运行是读取的共享内存中，在之后计算旋转矩阵时直接读取共享内存中的数据相比从
        全局内存（Global Memory）中读取速度大大提高。所以计算旋转矩阵任务划分的依据就在于每个$block$共享
        内存的大小。
        通过旋转矩阵对体数据进行旋转的操作与上述过程类似，不过可以不用将将旋转矩阵的内容读取到共享内存中
        ，因为旋转过程中这些数据只使用了一次，并不需要多次读取，所以这个操作也没有意义。
        \begin{figure}[htpb]
            \centering
            \includegraphics[width=3in]{5-7.jpg} %文件名中不能出现括号
            \caption{双线性插值需要保存的四个体素}
            \label{fig:5-7}
        \end{figure}
        \begin{figure}[htpb]
            \centering
            \includegraphics[width=3in]{5-8.jpg} %文件名中不能出现括号
            \caption{计算旋转矩阵任务划分示意图}
            \label{fig:5-8}
        \end{figure}

        正向投影时投影矩阵每一行中存储一条射线相关的体素，所以很自然的会想到一行数据对应到一个$block$，每
        个元素对应于一个$thread$，这样的任务划分会导致两个问题：一个是由于每个$thread$操作中所需的寄存器
        个数较多，如果同一个$block$中的$thread$个数太多时实际同时可运行的$thread$个数有限，第二个问题在于
        由于每个射线所经过的体素的个数并不相同，实际存储时投影矩阵中每一行都会有一部分空闲的存储区域，某
        些$thread$用于对这些空闲区域进行操作，这些操作的所需指令与普通$thread$的指令完全不同，而在CUDA中
        ，只有一个$warp$中的所有$thread$的指令相同才能达到最大的并行效率，否则有几组不同指令的$thread$便
        需要几次等待SM发生指令，从而不能达到最大的并行效率。所以在并行的时候尽量把不同操作的$thread$分类
        且分配到不同的$block$中，但实际中不能分别为不同的$block$分别设置大小，所以我们可以讲一行中的数据
        分成若干个$block$，这样至少可以保证前面几个$block$中的$thread$具有相同的指令发生过程，从而达到更
        高的并行效率，如图\ref{fig:5-9}所示，当不划分时一个$block$中包含不同类型的$thread$的概率为1，划分
        之后此概率降低为1/5。如果$block$分的越细$block$中$thread$一致的概率越高，但是也意味着$block$所绑定的SM能
        够同时运行的$thread$越少，也会影响并行的效率，所以此处$block$个数的划分需要一定的权衡，如表
        \ref{tb:5-3}所示为利用CUDA对单层重建进行加速后的$block$中包含不同数目$thread$时重建的时间，经过试
        验确定当一个$block$中的$thread$个数为128时并行效率最高，当$thread$个数小于128时由于$block$内
        $thread$个数不够无法充分并行运算，$thread$个数大于128时，由于某些$warp$内的$thread$指令不同也会影
        响并行速度，但是没有$block$内$thread$数目不足时的影响明显。
        \begin{table}[htpb]
            \centering
            \caption{$block$不同维度下CUDA重建时间（单位：ms）}
            \begin{tabular}{|l|*{4}{c|}}\hline
                \backslashbox{$block$维数}{投影个数} & 30 & 60 & 90 \\
                \hline
                32&174.435&341.230&512.975\\\hline
                64&169.888&332.131&500.840\\\hline
                128&164.638&324.124&486.454\\\hline
                256&166.319&322.980&490.014\\ \hline
            \end{tabular}
            \label{tb:5-3}
        \end{table}
        \begin{figure}[htpb]
            \centering
            \subfigure[不分割block]
            {
                \label{fig:5-9subfig:1}
                \includegraphics[width=3.5in]{5-9-1.jpg} %文件名中不能出现括号
            }
            \subfigure[分割block]
            {
                \label{fig:5-9subfig:2}
                \includegraphics[width=3.5in]{5-9-2.jpg} %文件名中不能出现括号
            }
            \caption{每一行矩阵元素任务划分示意图}
            \label{fig:5-9}
        \end{figure}

        反投影的过程是按照体素进行遍历，所以这个过程并行化CUDA代码的任务划分与计算投影矩阵类似，这里主要
        需要考虑的内容在于寄存器的个数均衡问题，主要的数据只需读取一遍，无需考虑使用共享内存对访存进行加速
        的优化，只需将循环遍历转换成并行算法即可得到较好加速效果。
        \subsection{实验结果}
        实验中采用的PC配置为Intel(R) Core(TM)2 Quad Q8400(2.66GHz)处理器，4G内存，NVIDIA公司的Geforce GTX
        560Ti显卡，该显卡的CUDA计算能力为2.1，具有48个SM，384个流处理器，显存容量为993MB，显存带宽为
        849.41MB/s 显存频率为2100.00MHz，SP核心频率为1.90GHz，开发环境为Visual Studio 2010 Ultimate加
        CUDA4.0开发环境。采用3D Shepp-Logan模型的投影数据进行实验。几何参数如图\ref{tb:5-4}所示，对锥束投
        影数据进行SART重建，比较分别在CPU和GPU加速情况下的重建结果。重建图像如图\ref{fig:5-10}所示。
        \begin{table}[htpb]
            \centering
            \caption{锥束扫描参数}
            \begin{tabular}{|l|*{2}{c|}}\hline
                扫描方式 & 圆周扫描  \\ \hline
                焦距(mm) & 50  \\ \hline
                探距(mm) & 1000  \\ \hline
                探测器长度(mm) & 128  \\ \hline
                探测器分辨率 & $256\times 256$  \\ \hline
                重建物体分辨率 & $256\times 256\times 256$  \\ \hline
                采集投影间隔($^o$) & 2 \\ \hline
                采集投影数(帧) & 180 \\ \hline
            \end{tabular}
            \label{tb:5-4}
        \end{table}
        \begin{figure}[htpb]
            \centering
            \subfigure[]
            {
                \label{fig:5-10subfig:1}
                \begin{minipage}{0.3\textwidth}
                    \includegraphics[width=2in]{5-10-1.jpg}\\
                    \vspace{0.01in}
                    \includegraphics[width=2in]{5-10-3.jpg}\\
                    \vspace{0.01in}
                    \includegraphics[width=2in]{5-10-5.jpg}
                \end{minipage}%
            }
            \hspace{0.5in}
            \subfigure[]
            {
                \label{fig:5-10subfig:2}
                \begin{minipage}{0.3\textwidth}
                    \includegraphics[width=2in]{5-10-2.jpg} \\
                    \vspace{0.01in}
                    \includegraphics[width=2in]{5-10-4.jpg} \\
                    \vspace{0.01in}
                    \includegraphics[width=2in]{5-10-6.jpg}
                \end{minipage}
            }
            \caption{z=100 y=100 x=125处(a)CPU(b)GPU迭代重建结果比较}
            \label{fig:5-10}
        \end{figure}
        为了评价原始SART重建的结果与经过GPU加速后的SART重建结果之间是否映入差异，我们还是用上文中提到的三
        种评价标准对两种方法的重建结果相似度进行比较。从表\ref{tb:5-5}可以看出利用GPU对SART算法进行加速时只
        是引入了极小的误差，主要时在任务分块时边界元素的插值元素不完备。两者重建的差别肉眼几乎不可分辨。
        \begin{table}[htpb]
            \centering
            \caption{原始SART重建与并行SART重建结果各项指标的比较}
            \begin{tabular}{|l|*{4}{c|}}\hline
                \backslashbox{评价指标}{投影个数} & z=100 & y=100&x=125 \\\hline
                重建误差d& 0.1220 & 0.1236 & 0.1443 \\
                \hline
                重建误差r& 0.0912& 0.0861 & 0.0824 \\
                \hline
                图像相似性系数$\varepsilon$& 0.9935 &0.9967&0.9928\\
                \hline
            \end{tabular}
            \label{tb:5-5}
        \end{table}

        在实验的不同阶段加速情况如表\ref{tb:5-6}所示，
        \begin{table}[htpb]
            \centering
            \caption{SART重建不同阶段所需时间对比表(单位：ms)}
            \begin{tabular}{|l|*{4}{c|}}\hline
            重建阶段& 原始算法 & GPU算法 & 加速比 \\\hline
            正向旋转&277.899&4.64016&59.89\\\hline
            正投影&559.732&30.8680&18.13\\\hline
            反投影&1124.31&29.0491&38.73\\\hline
            逆旋转&352.320&20.9081&16.85\\\hline
            \end{tabular}
            \label{tb:5-6}
        \end{table}

        在实际不同投影个数重建的时间与加速比如表\ref{tb:5-7}所示，
        \begin{table}[htpb]
            \centering
            \caption{各投影个数时SART串/并行重建时间对比}
            \begin{tabular}{|l|*{4}{c|}}\hline
            投影个数& 原始算法(秒) & GPU算法(秒)& 加速比 \\\hline
            30& 68.8180& 2.861&24\\\hline
            60& 138.803&5.396&25.72\\\hline
            90& 211.35& 7.953&26.57\\\hline
            \end{tabular}
            \label{tb:5-7}
        \end{table}

        从以上的实验结果表明利用GPU对SART算法做并行化计算，在引入极小误差的情况下重建的速度得到大大的提高
        ，基本上具有了在工程应用中的应用的可能。

        \section{本章小结}
        本章中首先对SART算法进行改进，利用投影系统刚体旋转的特性将投影系统的旋转映射成重建区域数据的旋转
        ，这样可以避免每个角度投影矩阵的计算，在引入极小误差的情况下能够提高重建的速度，并且使得SART算法
        更加适合于用CUDA并行化，在此基础上，我们发现此时的投影系统是关于射线中心轴对称的系统，所以在计算
        投影矩阵时可以只计算1/4大小射线对应的投影矩阵，其他3/4的投影矩阵可以通过对称性所得。进过如此处理
        之后的重建过程所需内存大大减小，从而使得在显存有限的图形卡上能够进行更高分别率的重建。之后就是对
        改进后的SART算法进行并行化改进，除了把原来程序中的循环改成对于各自的$thread$之外，主要的工作在于
        对各个阶段任务的划分，这需要熟悉CUDA中的GPU构架，设计出适合并行优化的任务划分方案以及访存策略。在
        完成以上工作之后通过实验证明经过GPU并行化的SART重建算法在引入极小的误差的情况下重建的速度得到极大
        的提高，可以满足实际工程的需求，达到了预期的目的。
    \chapter{总结与展望}
        \section{研究工作总结}
        锥束CT是临床介入手术中常用的设备，能够利用锥束投影数据直接进行体数据的重建，与传统的螺旋CT相比具
        有较高的垂直方向分辨率以及更高的数据利用率，而且采集过程中曝光时间较短，能够尽量降低辐射对于病人
        的影响，也保护介入医生的身体。锥束CT常见的重建方法分为解析重建和迭代重建两种方法，解析法利用基于
        Radon反变换的原理对三维区域进行重建，迭代法与之相比可以通过更少的投影个数达到相似的重建结果，但是
        迭代重建算法本身具有重建速度较慢的缺陷，然而随着GPU通用计算技术（GPGPU）的发展，可以将迭代算法实
        现并行化改造，从而实现重建速度的巨大提高成为一种可在实际工程中可用的重建算法。针对以上的背景，本
        文主要讨论了以下几方面的内容：

        (1)首先介绍了二维及三维重建的基本原理，分别从原理上研究了解析重建法与迭代重建法各自的特点以及重
        建各自所需的条件，最后介绍了重建结果的客观评价标准，为下面进一步的研究工作奠定理论基础。

        (2)对迭代重建中投影矩阵的计算进行研究，主要的研究内容包括不同的投影矩阵表示方法以及插值核的选择
        对重建结果的影响。之后对三种经典的迭代重建算法进行进一步的研究，最后通过各个实验的结果表明SART重
        建算法相对于其他两种经典的重建算法不仅对噪声投影图像重建的鲁棒性上而且重建的收敛速度上都具有一定的优势，所以
        在之后的研究中主要以SART重建算法为主要研究对象。

        (3)在锥束SART重建计算投影矩阵时，主要有两种计算的策略：基于体素驱动以及基于射线驱动，首先介绍了
        两种投影矩阵计算各自原理以及实现方法，基于体素驱动提出了一系列近似计算方法，理论上表明提出的这些
        近似对于重建的结果影响极小，而基于射线驱动的本身并不引入任何误差，且在较高分辨率重建时按照射线变
        遍历的射线驱动比体素驱动的方法速度更快。

        (4)利用不同投影的旋转是刚体变换的原理，将这个旋转映射成重建区域数据的旋转，基于这个原理对原始的
        SART进行一定的改进，不但在引入极小误差的同时可以提高重建的速度，而且使得改进后的算法更加适合于并
        行化加速。随后利用NVIDIA公司的CUDA并行开发环境对此算法进行并行化改进，极大地提高了SART重建的速度
        ，并且能够保证重建的质量。
        \section{展望}
        由于本人才疏学浅能力有限，在有限的硕士求学阶段中很多研究内容都是浅尝辄止，并没有进行很深入的研究
        。本课题还可以从以下的几个方面进行深入的挖掘：

        (1)对SART算法系列的改进算法如OS-SART等进行深入研究，在最新的研究中基于SART的改进算法具有各自不
        同的新的特性。

        (2)本文中没有对投影矩阵计算进行并行化算法的研究，而这其实是迭代重建非常重要的一个步骤，如何有效
        地对投影矩阵计算进行并行化是一个非常有意义的研究内容。

        (3)在现存的图形卡中，由于显存大小的问题，无法进行$256\times 256\times 256$更高分辨率的重建，所
        以对感兴趣区域（ROI）进行较高分辨率重建，而对其他部分进行较粗略的重建，这种多分辨率重建的方法也是
        一个很有意义的研究方向。

        (4)在CUDA对程序进行任务划分的时候仍然有很多优化的细节没有完全考虑，如何SART的CUDA重建算法进行更
        优化的任务划分也具有一定的研究价值。

        (5)在研究的过程中，对重建图像的后期处理也有一定研究，在可视化方面初步实现了重建图像的表面绘制和
        体绘制，但是并无创新之处，如何利用CUDA技术对重建图像可视化进行加速也是一项重要的功能。
\end{Main} % 结束正文
\begin{Thanks}
\end{Thanks}
\bibliographystyle{seuthesis}
\bibliography{seuthesis}
\end{document}
